fMRI Motivational Incentive Integration Brain Analyses
Debbie Yee
12/18/2020
- 1 Data Formatting and Scrubbing
- 1.1 Packages
- 1.2 Data Directories
- 1.3 Setting Subject IDs and Parameters
- 1.4 Read in Block and EV betas
- 1.5 Read in the behavioral & self-report data, format/average data based on relevant conditions
- 1.6 Summarising the incentive data based on relevant conditions
- 1.7 Summarising Bundled Beta Values in dorsal Anterior Cingulate Cortex (Knot 2 & 3, 4-6 ms)
- 1.8 Summarising Bundled Beta Values in Ventromedial PFC
- 1.9 Summarising Bundled Beta Values in Nucleus Accumbens
- 1.10 Summarising Bundled Beta Values in Caudate
- 1.11 Summarising Bundled Beta Values in Putamen
- 1.12 Summarising Bundled Beta Values in Anterior Insula
- 2 Motivational Incentive Integration Effects: Reward Rate and Self-Reported Motivation Ratings
- 2.1 Extended Data Figure 1-2: RT between baseline and incentive block
- 2.2 Extended Data Paired T-test: RT between baseline and incentive block
- 2.3 Extended Data Figure 1-2: Accuracy between baseline and incentive block
- 2.4 Extended Data Paired T-test: Accuracy between baseline and incentive block
- 2.5 Extended Data Figure 1-2: RT switch costs between baseline and incentive block
- 2.6 Extended Data Paired T-test: RT switch costs between baseline and incentive block
- 2.7 Plot RT by switch and non-switch trials
- 2.8 Extended Data Save Figure 1-2: Combine RT, Accuracy, and RT switch cost plots
- 2.9 Reward Rate by Money and Liquid Incentives
- 2.10 RT and Accuracy by Money and Liquid Incentives
- 2.11 Extended Data Table 1-2: RT and Accuracy by Money and Liquid Incentives
- 2.12 Extended Table 1-3: Multiple Regression Between Reward Rate and Motivation Ratings
- 2.13 Hierarchical Regression with Reward Rate and Motivation Ratings
- 2.14 Reward Rate by Motivation Ratings By Subject (Scatterplot)
- 2.15 Reward Rate by Motivation Ratings by Monetary Reward and Liquid (Scatterplot)
- 2.16 Motivation Ratings by Monetary Reward and Liquid (Barplot)
- 2.17 Reward Rate by Monetary Reward and Liquid (Barplot)
- 2.18 Manuscript Save Figure 1
- 2.19 Extended Data Figure 1-4: RT and Accuracy by Monetary Reward and Liquid (Barplot)
- 2.20 Extended Data Save Figure 1-4
- 3 Dorsal Anterior Cingulate Cortex (dACC) Encodes Both Monetary and Liquid Incentives in Bundled Beta Estimates
- 3.1 Block Beta Values by Liquid in dACC
- 3.2 Event-Related Beta Values by Money in dACC
- 3.3 Bundled Beta Values by Money in dACC
- 3.4 Table of the Block, Event-Related,and Bundled Models
- 3.5 Manuscript Figure 2: Plotting dACC Blocked (Liquid) Effects & Cue (Money) Effects
- 3.6 Manuscript Figure 2: Plotting dACC Bundled Effects
- 3.7 Manuscript Figure 2: Draw Tentzero function of trialwise beta estimates
- 3.8 Manuscript Save Figure 2
- 4 Reward Rate and Self-Report Motivation Ratings Predicted by Dorsal ACC Bundled Betas
- 5 Dorsal ACC Effects on Reward Rate Mediated by Self-Report Motivation Ratings (Within-Subjects Mediation Analyses)
- 5.1 Format the Data for Within-Subjects Mediation Analysis
- 5.2 Step 0: Motivation Conditions Predict Reward Rate (Y), dACC Betas, and Motivation Ratings
- 5.3 Step 1: Reward Rate (Y) by dACC Betas (X)
- 5.4 Step 2: Motivation Ratings (M) Predicted by dACC Betas
- 5.5 Step 3: dACC Betas Predict Reward Rate Controling for Motivation Ratings
- 5.6 Plot all 3 Models
- 5.7 Step 4: Run the mediation using bmlm package
- 5.8 Inverse mediation, with RR mediating effects of dACC on motivation ratings.
- 6 Between-Subjects correlations between Reward Rate, Motivation Ratings, and dACC
- 7 Dorsal ACC Selectively Encodes Subjective Motivational Value and Modulates Motivated Task Performance
- 8 Exploratory Visulization and Plots (Not included in publication)
- 9 ROI Analyses: Putamen
- 10 Exploratory Visualization: Congruency
- 11 Exploratory Visualization: Reward Rate switch costs
- 12 Exploratory Visualization: RT Switch Costs
1 Data Formatting and Scrubbing
1.1 Packages
1.2 Data Directories
1.3 Setting Subject IDs and Parameters
1.4 Read in Block and EV betas
# Schaefer Parcels (400)
betas.block.sch<-read.csv(file = paste0(data.path.betas,"betas_block_schaefer_long.csv")) %>%
dplyr::select(-X) %>%
mutate(Liquid=factor(block, levels = c("BASE","JUICE","NEUT","SALT"),
labels = c("Baseline","Juice","Neutral","Saltwater")))
betas.ev.sch<-read.csv(file = paste0(data.path.betas,"betas_ev_schaefer_long.csv")) %>%
dplyr::select(-X) %>%
separate(col = condition, into = c("Liquid","Reward"), sep = "_", remove = FALSE) %>%
mutate(tentplot = tent+1, AllvsNoReward=ifelse(Liquid=="Noreward", yes = "NoReward", no = "Reward"))
# Subcortical Parcels (19)
betas.block.sub<-read.csv(file = paste0(data.path.betas,"betas_block_subcortical_long.csv")) %>%
dplyr::select(-X) %>%
mutate(Liquid=factor(block, levels = c("BASE","JUICE","NEUT","SALT"),
labels = c("Baseline","Juice","Neutral","Saltwater")))
betas.ev.sub<-read.csv(file = paste0(data.path.betas,"betas_ev_subcortical_long.csv")) %>%
dplyr::select(-X) %>%
separate(col = condition, into = c("Liquid","Reward"), sep = "_", remove = FALSE) %>%
mutate(tentplot = tent+1, AllvsNoReward=ifelse(Liquid=="Noreward", yes = "NoReward", no = "Reward"))1.5 Read in the behavioral & self-report data, format/average data based on relevant conditions
#data.beta.block<-read.csv(paste0(data.path.betas,"betas_block.csv")) %>% select(-X)
#data.beta.EV<-read.csv(paste0(data.path.betas,"betas_EV.csv")) %>% select(-X)
data.behavior <- read.delim(data.path.behavior)
data.runkey <- read.csv(runkey.path)
data.selfreport<-read.csv(selfreport.path)
# BASELINE subject data (filter using data pipeline from dplyr)
base <- filter(data.behavior,runType=="baseline") %>%
dplyr::select(subID,runID,runType,trial,RT,ACC,switch_trial,task,congruency) %>%
mutate(ERR=ifelse(ACC==0,1,0))
base.means <- base %>% group_by(subID) %>%
summarise(n=n(), meanACC = mean(ACC, na.rm=TRUE), seACC = sd(ACC, na.rm=TRUE)/sqrt(n-1),
meanERR = mean(ERR, na.rm=TRUE), seERR = sd(ERR,na.rm=TRUE)/sqrt(n-1))
baseRT.means <- base %>% group_by(subID) %>% filter(ACC==1) %>%
summarise(n=n(), meanRT = mean(RT), seRT = sd(RT)/sqrt(n-1))
# INCENTIVE subject data (filter using data pipeline from dplyr)
incentive <- data.behavior %>% filter(runType=="incentive" & !is.na(RT)) %>%
dplyr::select(subID,runID,runType,trial,RT,ACC,switch_trial,task,congruency,Feedback,
moneyreward,money = moneysymbol) %>%
inner_join(data.runkey, by=c("subID","runID")) %>%
mutate(critRT = critRT,
subRewarded = (RT<critRT & ACC==1)*1,
ERR=ifelse(ACC==0,1,0),
moneyCode = factor(moneyreward, levels = c(1,2,4), labels = c(-1,0,1)),
liqCode = factor(liquid,levels=c("saltwater","neutral","juice"),
labels=c(-1,0,1)),
rankCode = factor(rank,levels=c(1,2,3),labels=c(1,0,-1)),
orderCode = factor(rank,levels=c(1,2,3),labels=c(-1,0,1)),
liqCodeJvN = factor(liquid,levels=c("saltwater","neutral","juice"),
labels=c(0,-1,1)),
liqCodeSvN = factor(liquid,levels=c("saltwater","neutral","juice"),
labels=c(-1,1,0)),
subID=as.factor(subID)) %>%
group_by(subID) %>%
mutate(zRT = scale(RT,center = TRUE, scale = TRUE))
incentive$subID<-as.numeric(levels(incentive$subID))[incentive$subID]
incentive$moneyCode<-as.numeric(levels(incentive$moneyCode)[incentive$moneyCode])
incentive$liqCode<-as.numeric(levels(incentive$liqCode)[incentive$liqCode])
incentive$rankCode<-as.numeric(levels(incentive$rankCode)[incentive$rankCode])
incentive$orderCode<-as.numeric(levels(incentive$orderCode)[incentive$orderCode])
incentive$liqCodeJvN<-as.numeric(levels(incentive$liqCodeJvN)[incentive$liqCodeJvN])
incentive$liqCodeSvN<-as.numeric(levels(incentive$liqCodeSvN)[incentive$liqCodeSvN])1.6 Summarising the incentive data based on relevant conditions
# summarise means of incentive data for each subject, across all conditions
incentive.means <- incentive %>% filter(!subID %in% sub.exclude) %>% group_by(subID) %>%
summarise(n=n(), meanRR = mean(subRewarded, na.rm=TRUE),
seRR = sd(subRewarded, na.rm=TRUE)/sqrt(n-1),
meanACC = mean(ACC, na.rm=TRUE), seACC = sd(ACC, na.rm=TRUE)/sqrt(n-1),
meanERR = mean(ERR, na.rm=TRUE), seERR = sd(ERR,na.rm=TRUE)/sqrt(n-1))
incentiveRT.means <- incentive %>% group_by(subID) %>% filter(!subID %in% sub.exclude,ACC==1) %>%
summarise(n=n(), meanRT = mean(RT, na.rm=TRUE), seRT = sd(RT)/sqrt(n-1))
# summarise by subject and congruency across all conditions
incentive.means.cong <- incentive %>% filter(!subID %in% sub.exclude) %>% group_by(subID,congruency) %>%
summarise(n=n(), meanRR = mean(subRewarded, na.rm=TRUE),
seRR = sd(subRewarded, na.rm=TRUE)/sqrt(n-1),
meanACC = mean(ACC, na.rm=TRUE), seACC = sd(ACC, na.rm=TRUE)/sqrt(n-1),
meanERR = mean(ERR, na.rm=TRUE), seERR = sd(ERR,na.rm=TRUE)/sqrt(n-1))
incentiveRT.means.cong <- incentive %>% group_by(subID,congruency) %>% filter(!subID %in% sub.exclude,ACC==1) %>%
summarise(n=n(), meanRT = mean(RT, na.rm=TRUE), seRT = sd(RT)/sqrt(n-1))
# summarise by subject and switch cost across all conditions
incentive.means.switch <- incentive %>% filter(!subID %in% sub.exclude, trial!=1) %>% group_by(subID,switch_trial) %>%
summarise(n=n(), meanRR = mean(subRewarded, na.rm=TRUE),
seRR = sd(subRewarded, na.rm=TRUE)/sqrt(n-1),
meanACC = mean(ACC, na.rm=TRUE), seACC = sd(ACC, na.rm=TRUE)/sqrt(n-1),
meanERR = mean(ERR, na.rm=TRUE), seERR = sd(ERR,na.rm=TRUE)/sqrt(n-1))
incentiveRT.means.switch <- incentive %>% group_by(subID,switch_trial) %>% filter(!subID %in% sub.exclude,ACC==1,trial!=1) %>%
summarise(n=n(), meanRT = mean(RT, na.rm=TRUE), seRT = sd(RT)/sqrt(n-1))
# summarise the means of the incentive data for each subject, grouped by condition
incentive9RT.means <- incentive %>% filter(!subID %in% sub.exclude) %>% group_by(subID, liquid, money) %>%
filter(!is.na(subRewarded), ACC==1) %>% summarise(n=n(), meanRT = mean(RT)) %>%
ungroup(liquid) %>% dplyr::rename(Liquid=liquid) %>%
mutate(Liquid = factor(Liquid, levels = c("juice","neutral","saltwater"), labels = c("Juice","Neutral","Saltwater")),
Reward = as.character(factor(money,levels=c("$","$$","$$$$"),labels=c(1,2,4))),
MoneyLabel = factor(x = money, levels = c("$","$$","$$$$"), labels = c("Low", "Med", "High")),
moneyCode = as.numeric(as.character(factor(money,levels=c("$","$$","$$$$"),labels=c(-1,0,1)))),
liqCode = as.numeric(as.character(factor(Liquid,levels=c("Juice","Neutral","Saltwater"),
labels=c(1,0,-1))))) %>%
ungroup(subID)
incentive9.means = incentive %>% filter(!subID %in% sub.exclude, !is.na(subRewarded)) %>% group_by(subID, liquid, money) %>%
summarise(n = n(), meanRR = mean(subRewarded), meanACC = mean(ACC), meanERR = mean(ERR)) %>%
ungroup(liquid) %>% dplyr::rename(Liquid=liquid) %>%
mutate(Liquid = factor(Liquid, levels = c("juice","neutral","saltwater"), labels = c("Juice","Neutral","Saltwater")),
Reward = as.character(factor(money,levels=c("$","$$","$$$$"),labels=c(1,2,4))),
MoneyLabel = factor(x = money, levels = c("$","$$","$$$$"), labels = c("Low", "Med", "High")),
moneyCode = as.numeric(as.character(factor(money,levels=c("$","$$","$$$$"),labels=c(-1,0,1)))),
liqCode = as.numeric(as.character(factor(Liquid,levels=c("Juice","Neutral","Saltwater"),labels=c(1,0,-1))))) %>%
ungroup(subID) %>%
left_join(incentive9RT.means, by = c("subID","money","Liquid","Reward","MoneyLabel","moneyCode","liqCode"))
# summarise the means of the incentive data for each subject, grouped by congruency and condition
incentive9RT.means.cong <- incentive %>% filter(!subID %in% sub.exclude) %>% group_by(subID, congruency,liquid, money) %>%
filter(!is.na(subRewarded), ACC==1) %>% summarise(n=n(), meanRT = mean(RT)) %>%
ungroup(liquid) %>% dplyr::rename(Liquid=liquid) %>%
mutate(Liquid = factor(Liquid, levels = c("juice","neutral","saltwater"), labels = c("Juice","Neutral","Saltwater")),
Reward = as.character(factor(money,levels=c("$","$$","$$$$"),labels=c(1,2,4))),
MoneyLabel = factor(x = money, levels = c("$","$$","$$$$"), labels = c("Low", "Med", "High")),
moneycode = as.numeric(as.character(factor(money,levels=c("$","$$","$$$$"),labels=c(-1,0,1)))),
liqcode = as.numeric(as.character(factor(Liquid,levels=c("juice","neutral","saltwater"),
labels=c(1,0,-1))))) %>%
ungroup(subID)
incentive9.means.cong = incentive %>% filter(!subID %in% sub.exclude, !is.na(subRewarded)) %>%
group_by(subID, congruency, liquid, money) %>%
summarise(n = n(), meanRR = mean(subRewarded), meanACC = mean(ACC), meanERR = mean(ERR)) %>%
ungroup(liquid) %>% dplyr::rename(Liquid=liquid) %>%
mutate(Liquid = factor(Liquid, levels = c("juice","neutral","saltwater"), labels = c("Juice","Neutral","Saltwater")),
Reward = as.character(factor(money,levels=c("$","$$","$$$$"),labels=c(1,2,4))),
MoneyLabel = factor(x = money, levels = c("$","$$","$$$$"), labels = c("Low", "Med", "High")),
moneyCode = as.numeric(as.character(factor(money,levels=c("$","$$","$$$$"),labels=c(-1,0,1)))),
liqCode = as.numeric(as.character(factor(Liquid,levels=c("Juice","Neutral","Saltwater"),labels=c(1,0,-1))))) %>%
ungroup(subID)
# summarise the means of the incentive data for each subject, grouped by switch cost and condition
incentive9RT.means.switch <- incentive %>% filter(!subID %in% sub.exclude) %>%
group_by(subID, switch_trial,liquid,money) %>%
filter(!is.na(subRewarded), ACC==1, trial!=1) %>% summarise(n=n(), meanRT = mean(RT)) %>%
ungroup(liquid) %>% dplyr::rename(Liquid=liquid) %>%
mutate(Liquid = factor(Liquid, levels = c("juice","neutral","saltwater"), labels = c("Juice","Neutral","Saltwater")),
Reward = as.character(factor(money,levels=c("$","$$","$$$$"),labels=c(1,2,4))),
MoneyLabel = factor(x = money, levels = c("$","$$","$$$$"), labels = c("Low", "Med", "High")),
moneycode = as.numeric(as.character(factor(money,levels=c("$","$$","$$$$"),labels=c(-1,0,1)))),
liqcode = as.numeric(as.character(factor(Liquid,levels=c("Juice","Neutral","Saltwater"),
labels=c(1,0,-1)))),
switchcode = as.numeric(as.character(factor(switch_trial, levels = c("switch","repeat"), labels = c(1,0))))) %>%
ungroup(subID)
incentive9.means.switch = incentive %>% filter(!subID %in% sub.exclude, !is.na(subRewarded), trial!=1) %>%
group_by(subID, switch_trial, liquid, money) %>%
summarise(n = n(), meanRR = mean(subRewarded), meanACC = mean(ACC), meanERR = mean(ERR)) %>%
ungroup(liquid) %>% dplyr::rename(Liquid=liquid) %>%
mutate(Liquid = factor(Liquid, levels = c("juice","neutral","saltwater"), labels = c("Juice","Neutral","Saltwater")),
Reward = as.character(factor(money,levels=c("$","$$","$$$$"),labels=c(1,2,4))),
MoneyLabel = factor(x = money, levels = c("$","$$","$$$$"), labels = c("Low", "Med", "High")),
moneyCode = as.numeric(as.character(factor(money,levels=c("$","$$","$$$$"),labels=c(-1,0,1)))),
liqCode = as.numeric(as.character(factor(Liquid,levels=c("Juice","Neutral","Saltwater"),labels=c(1,0,-1))))) %>%
ungroup(subID)
# summarise the self report data
# motivation ratings
selfreport.motivation <- data.selfreport %>% filter(!subID %in% sub.exclude) %>%
dplyr::select(subID, juice1_motivetrial,juice2_motivetrial,juice4_motivetrial,
neut1_motivetrial,neut2_motivetrial,neut4_motivetrial,
salt1_motivetrial,salt2_motivetrial,salt4_motivetrial) %>%
gather(key = condition, value = motive_rating, juice1_motivetrial:salt4_motivetrial) %>%
separate(col = condition, into = c("condition","misc"), sep="_", remove = TRUE) %>%
dplyr::select(-misc) %>%
mutate(condition = replace(condition, list = which(condition=="juice1"), values = "Juice_$"),
condition = replace(condition, list = which(condition=="juice2"), values = "Juice_$$"),
condition = replace(condition, list = which(condition=="juice4"), values = "Juice_$$$$"),
condition = replace(condition, list = which(condition=="neut1"), values = "Neutral_$"),
condition = replace(condition, list = which(condition=="neut2"), values = "Neutral_$$"),
condition = replace(condition, list = which(condition=="neut4"), values = "Neutral_$$$$"),
condition = replace(condition, list = which(condition=="salt1"), values = "Saltwater_$"),
condition = replace(condition, list = which(condition=="salt2"), values = "Saltwater_$$"),
condition = replace(condition, list = which(condition=="salt4"), values = "Saltwater_$$$$")) %>%
separate(col = condition, into = c("Liquid","money"), sep = "_") %>%
mutate(Reward = as.character(factor(money,levels=c("$","$$","$$$$"),labels=c(1,2,4))))
# liking ratings
selfreport.liking <- data.selfreport %>% filter(!subID %in% sub.exclude) %>%
dplyr::select(subID, juice1_liketrial,juice2_liketrial,juice4_liketrial,
neut1_liketrial,neut2_liketrial,neut4_liketrial,
salt1_liketrial,salt2_liketrial,salt4_liketrial) %>%
gather(key = condition, value = liking_rating, juice1_liketrial:salt4_liketrial) %>%
separate(col = condition, into = c("condition","misc"), sep="_", remove = TRUE) %>%
dplyr::select(-misc) %>%
mutate(condition = replace(condition, list = which(condition=="juice1"), values = "Juice_$"),
condition = replace(condition, list = which(condition=="juice2"), values = "Juice_$$"),
condition = replace(condition, list = which(condition=="juice4"), values = "Juice_$$$$"),
condition = replace(condition, list = which(condition=="neut1"), values = "Neutral_$"),
condition = replace(condition, list = which(condition=="neut2"), values = "Neutral_$$"),
condition = replace(condition, list = which(condition=="neut4"), values = "Neutral_$$$$"),
condition = replace(condition, list = which(condition=="salt1"), values = "Saltwater_$"),
condition = replace(condition, list = which(condition=="salt2"), values = "Saltwater_$$"),
condition = replace(condition, list = which(condition=="salt4"), values = "Saltwater_$$$$")) %>%
separate(col = condition, into = c("Liquid","money"), sep = "_") %>%
mutate(Reward = as.character(factor(money,levels=c("$","$$","$$$$"),labels=c(1,2,4))))
# performance ratings
selfreport.perform <- data.selfreport %>% filter(!subID %in% sub.exclude) %>%
dplyr::select(subID, juice1_performtrial,juice2_performtrial,juice4_performtrial,
neut1_performtrial,neut2_performtrial,neut4_performtrial,
salt1_performtrial,salt2_performtrial,salt4_performtrial) %>%
gather(key = condition, value = perform_rating, juice1_performtrial:salt4_performtrial) %>%
separate(col = condition, into = c("condition","misc"), sep="_", remove = TRUE) %>%
dplyr::select(-misc) %>%
mutate(condition = replace(condition, list = which(condition=="juice1"), values = "Juice_$"),
condition = replace(condition, list = which(condition=="juice2"), values = "Juice_$$"),
condition = replace(condition, list = which(condition=="juice4"), values = "Juice_$$$$"),
condition = replace(condition, list = which(condition=="neut1"), values = "Neutral_$"),
condition = replace(condition, list = which(condition=="neut2"), values = "Neutral_$$"),
condition = replace(condition, list = which(condition=="neut4"), values = "Neutral_$$$$"),
condition = replace(condition, list = which(condition=="salt1"), values = "Saltwater_$"),
condition = replace(condition, list = which(condition=="salt2"), values = "Saltwater_$$"),
condition = replace(condition, list = which(condition=="salt4"), values = "Saltwater_$$$$")) %>%
separate(col = condition, into = c("Liquid","money"), sep = "_") %>%
mutate(Reward = as.character(factor(money,levels=c("$","$$","$$$$"),labels=c(1,2,4))))1.7 Summarising Bundled Beta Values in dorsal Anterior Cingulate Cortex (Knot 2 & 3, 4-6 ms)
#'dACC' parcel ids from the incentive integration analyses
parcels.dACC<-c(107,110,108,311,312)
# Looping over 'dACC' parcels and create data frame of combining effects from both parcels
for (pid in parcels.dACC) { # pid<-108
# subset data for each parcel
# blocked/sustained effects data
tmp.block <- betas.block.sch %>% filter(Liquid!="Baseline", parcel.id==pid) %>%
dplyr::rename(value.block = value) %>%
mutate(Liquid = as.character(Liquid))
# bundled effects & behavioral data for cue related activity
tmp <- betas.ev.sch %>% filter(condition!="Noreward", parcel.id==pid,
(tentplot>=TR_cue[1] & tentplot<=TR_cue[2])) %>%
group_by(subID,Liquid,Reward) %>%
dplyr::summarise(mean_value = mean(value)) %>% dplyr::rename(mean_value.ev = mean_value) %>%
mutate(liqCode = as.numeric(as.character(factor(Liquid,levels=c("Saltwater","Neutral","Juice"),labels=c(-1,0,1)))),
moneyCode = as.numeric(as.character(factor(Reward,levels=c("1","2","4"),labels=c(-1,0,1))))) %>%
left_join(y = tmp.block, by = c("subID","Liquid")) %>%
mutate(value.bundle = mean_value.ev + value.block) %>%
inner_join(y = incentive9.means, by = c("subID","Liquid","Reward","liqCode","moneyCode")) %>%
inner_join(y = selfreport.motivation, by = c("subID","Liquid","Reward","money")) %>%
inner_join(y = selfreport.liking, by = c("subID","Liquid","Reward","money")) %>%
inner_join(y = selfreport.perform, by = c("subID","Liquid","Reward","money"))
#tmp$liqCode <- as.numeric(levels(tmp$liqCode)[tmp$liqCode])
#tmp$moneyCode <- as.numeric(levels(tmp$moneyCode)[tmp$moneyCode])
if (pid == parcels.dACC[1]) {
data.dACC<-tmp
} else {
data.dACC<- rbind(data.dACC,tmp)
}
}`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)
`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)
`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)
`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)
`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)# Summarise the bundled values by Reward Rate across the parcels 107, 108, 110, 311, 312
data.dACC.sum <- data.dACC %>% group_by(subID, Liquid, Reward) %>%
summarise(value.block = mean(value.block),
value.ev = mean(mean_value.ev),
value.bundle.mean = mean(value.bundle),
meanRR = mean(meanRR),
meanERR = mean(meanRR),
meanACC = mean(meanACC),
meanRT = mean(meanRT), meanRTlog = log(meanRT),
motive_rating = mean(motive_rating), liking_rating = mean(liking_rating),
liqCode = mean(liqCode), moneyCode = mean(moneyCode)) %>%
ungroup(Liquid,Reward) %>% group_by(subID) %>%
mutate(MoneyLabel = factor(x = Reward, levels = c(1,2,4), labels = c("Low", "Med", "High")),
mot9Code = liqCode + moneyCode,
meanRR.normed = scale(meanRR, center = TRUE, scale = TRUE),
meanRR.centered = scale(meanRR, center = TRUE, scale = FALSE),
motive_rating.normed = scale(motive_rating, center = TRUE, scale = TRUE),
motive_rating.centered = scale(motive_rating, center = TRUE, scale = FALSE),
value.bundle.mean.normed = scale(value.bundle.mean, center = TRUE, scale = TRUE),
value.bundle.mean.centered = scale(value.bundle.mean, center = TRUE, scale = TRUE),
value.bundle.mean.mot9 = value.bundle.mean*mot9Code,
meanRR.mot9 = meanRR*mot9Code,
motive_rating.mot9 = motive_rating*mot9Code) %>%
unite(col = "Incentive9", Liquid:Reward, sep="_", remove = FALSE)`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)data.dACC.sum$motive_rating.normed[is.nan(data.dACC.sum$motive_rating.normed)]<-0
# Summarise by subject
data.dACC.sum.sub <- data.dACC %>% group_by(subID) %>%
summarise(value.block = mean(value.block),
value.ev = mean(mean_value.ev),
value.bundle.mean = mean(value.bundle),
meanRR = mean(meanRR),
meanERR = mean(meanRR), meanACC = mean(meanACC),
motive_rating = mean(motive_rating), liking_rating = mean(liking_rating),
liqCode = mean(liqCode), moneyCode = mean(moneyCode)) %>%
#ungroup(Liquid,Reward) %>% group_by(subID) %>%
mutate(meanRR.normed = scale(meanRR, center = TRUE, scale = TRUE),
motive_rating.normed = scale(motive_rating, center = TRUE, scale = TRUE),
value.bundle.mean.normed = scale(value.bundle.mean, center = TRUE, scale = TRUE))`summarise()` ungrouping output (override with `.groups` argument)1.8 Summarising Bundled Beta Values in Ventromedial PFC
#'vmpfc' parcel ids from the incentive integration analyses
#parcels.vmpfc<-c(169,379)
parcels.vmpfc<-c(168, 169, 173, 174, 379, 380, 381, 382)
# Looping over 'vmpfc' parcels and create data frame of combining effects from both parcels
for (pid in parcels.vmpfc) { # pid<-169
# subset data for each parcel
# blocked/sustained effects data
tmp.block <- betas.block.sch %>% filter(Liquid!="Baseline", parcel.id==pid) %>%
dplyr::rename(value.block = value) %>%
mutate(Liquid = as.character(Liquid))
# bundled effects & behavioral data for cue related activity
tmp <- betas.ev.sch %>% filter(condition!="Noreward", parcel.id==pid,
(tentplot>=TR_cue[1] & tentplot<=TR_cue[2])) %>%
group_by(subID,Liquid,Reward) %>%
summarise(mean_value = mean(value)) %>% dplyr::rename(mean_value.ev = mean_value) %>%
mutate(liqCode = as.numeric(as.character(factor(Liquid,levels=c("Saltwater","Neutral","Juice"),labels=c(-1,0,1)))),
moneyCode = as.numeric(as.character(factor(Reward,levels=c("1","2","4"),labels=c(-1,0,1))))) %>%
left_join(y = tmp.block, by = c("subID","Liquid")) %>%
mutate(value.bundle = mean_value.ev + value.block) %>%
inner_join(y = incentive9.means, by = c("subID","Liquid","Reward","liqCode","moneyCode")) %>%
inner_join(y = selfreport.motivation, by = c("subID","Liquid","Reward","money")) %>%
inner_join(y = selfreport.liking, by = c("subID","Liquid","Reward","money")) %>%
inner_join(y = selfreport.perform, by = c("subID","Liquid","Reward","money"))
#tmp$liqCode <- as.numeric(levels(tmp$liqCode)[tmp$liqCode])
#tmp$moneyCode <- as.numeric(levels(tmp$moneyCode)[tmp$moneyCode])
if (pid == parcels.vmpfc[1]) {
data.vmpfc<-tmp
} else {
data.vmpfc<- rbind(data.vmpfc,tmp)
}
}`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)
`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)
`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)
`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)
`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)
`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)
`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)
`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)# Summarise the bundled values across the vmpfc parcels 169 and 379
data.vmpfc.sum <- data.vmpfc %>% group_by(subID, Liquid, Reward) %>%
summarise(value.bundle.mean = mean(value.bundle), meanRR = mean(meanRR),
meanERR = mean(meanRR), meanACC = mean(meanACC),
motive_rating = mean(motive_rating), liking_rating = mean(liking_rating),
liqCode = mean(liqCode), moneyCode = mean(moneyCode)) %>%
ungroup(Liquid,Reward) %>% group_by(subID) %>%
mutate(mot9Code = liqCode + moneyCode,
meanRR.normed = scale(meanRR, center = TRUE, scale = TRUE),
meanRR.centered = scale(meanRR, center = TRUE, scale = FALSE),
motive_rating.normed = scale(motive_rating, center = TRUE, scale = TRUE),
motive_rating.centered = scale(motive_rating, center = TRUE, scale = FALSE),
value.bundle.mean.normed = scale(value.bundle.mean, center = TRUE, scale = TRUE),
value.bundle.mean.centered = scale(value.bundle.mean, center = TRUE, scale = TRUE),
value.bundle.mean.mot9 = value.bundle.mean*mot9Code,
meanRR.mot9 = meanRR*mot9Code,
motive_rating.mot9 = motive_rating*mot9Code)`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)data.vmpfc.sum$motive_rating.normed[is.nan(data.vmpfc.sum$motive_rating.normed)]<-0
# Summarise by subject
data.vmpfc.sum.sub <- data.vmpfc %>% group_by(subID) %>%
summarise(value.bundle.mean = mean(value.bundle), meanRR = mean(meanRR),
meanERR = mean(meanRR), meanACC = mean(meanACC),
motive_rating = mean(motive_rating), liking_rating = mean(liking_rating),
liqCode = mean(liqCode), moneyCode = mean(moneyCode)) %>%
#ungroup(Liquid,Reward) %>% group_by(subID) %>%
mutate(meanRR.normed = scale(meanRR, center = TRUE, scale = TRUE),
motive_rating.normed = scale(motive_rating, center = TRUE, scale = TRUE),
value.bundle.mean.normed = scale(value.bundle.mean, center = TRUE, scale = TRUE))`summarise()` ungrouping output (override with `.groups` argument)1.9 Summarising Bundled Beta Values in Nucleus Accumbens
# 'accumbens' parcel ids from the incentive integration analyses
parcels.accumbens<-c(350,351)
# Looping over 'vmpfc' parcels and create data frame of combining effects from both parcels
for (pid in parcels.accumbens) { # pid<-350
# subset data for each parcel
# blocked/sustained effects data
tmp.block <- betas.block.sub %>% filter(Liquid!="Baseline", parcel.id==pid) %>%
dplyr::rename(value.block = value) %>%
mutate(Liquid = as.character(Liquid))
# bundled effects & behavioral data for cue related activity
tmp <- betas.ev.sub %>% filter(condition!="Noreward", parcel.id==pid,
(tentplot>=TR_cue[1] & tentplot<=TR_cue[2])) %>%
group_by(subID,Liquid,Reward) %>%
summarise(mean_value = mean(value)) %>% dplyr::rename(mean_value.ev = mean_value) %>%
mutate(liqCode = as.numeric(as.character(factor(Liquid,levels=c("Saltwater","Neutral","Juice"),labels=c(-1,0,1)))),
moneyCode = as.numeric(as.character(factor(Reward,levels=c("1","2","4"),labels=c(-1,0,1))))) %>%
left_join(y = tmp.block, by = c("subID","Liquid")) %>%
mutate(value.bundle = mean_value.ev + value.block) %>%
inner_join(y = incentive9.means, by = c("subID","Liquid","Reward","liqCode","moneyCode")) %>%
inner_join(y = selfreport.motivation, by = c("subID","Liquid","Reward","money")) %>%
inner_join(y = selfreport.liking, by = c("subID","Liquid","Reward","money")) %>%
inner_join(y = selfreport.perform, by = c("subID","Liquid","Reward","money"))
#tmp$liqCode <- as.numeric(levels(tmp$liqCode)[tmp$liqCode])
#tmp$moneyCode <- as.numeric(levels(tmp$moneyCode)[tmp$moneyCode])
if (pid == parcels.accumbens[1]) {
data.accumbens<-tmp
} else {
data.accumbens<- rbind(data.accumbens,tmp)
}
}`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)
`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)# Summarise the bundled values across the accumbens parcels 350 and 351
data.accumbens.sum <- data.accumbens %>% group_by(subID, Liquid, Reward) %>%
summarise(value.bundle.mean = mean(value.bundle), meanRR = mean(meanRR),
meanERR = mean(meanRR), meanACC = mean(meanACC),
motive_rating = mean(motive_rating), liking_rating = mean(liking_rating),
liqCode = mean(liqCode), moneyCode = mean(moneyCode)) %>%
ungroup(Liquid,Reward) %>% group_by(subID) %>%
mutate(mot9Code = liqCode + moneyCode,
meanRR.normed = scale(meanRR, center = TRUE, scale = TRUE),
meanRR.centered = scale(meanRR, center = TRUE, scale = FALSE),
motive_rating.normed = scale(motive_rating, center = TRUE, scale = TRUE),
motive_rating.centered = scale(motive_rating, center = TRUE, scale = FALSE),
value.bundle.mean.normed = scale(value.bundle.mean, center = TRUE, scale = TRUE),
value.bundle.mean.centered = scale(value.bundle.mean, center = TRUE, scale = TRUE),
value.bundle.mean.mot9 = value.bundle.mean*mot9Code,
meanRR.mot9 = meanRR*mot9Code,
motive_rating.mot9 = motive_rating*mot9Code)`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)data.accumbens.sum$motive_rating.normed[is.nan(data.accumbens.sum$motive_rating.normed)]<-0
# Summarise by subject
data.accumbens.sum.sub <- data.accumbens %>% group_by(subID) %>%
summarise(value.bundle.mean = mean(value.bundle), meanRR = mean(meanRR),
meanERR = mean(meanRR), meanACC = mean(meanACC),
motive_rating = mean(motive_rating), liking_rating = mean(liking_rating),
liqCode = mean(liqCode), moneyCode = mean(moneyCode)) %>%
#ungroup(Liquid,Reward) %>% group_by(subID) %>%
mutate(meanRR.normed = scale(meanRR, center = TRUE, scale = TRUE),
motive_rating.normed = scale(motive_rating, center = TRUE, scale = TRUE),
value.bundle.mean.normed = scale(value.bundle.mean, center = TRUE, scale = TRUE))`summarise()` ungrouping output (override with `.groups` argument)1.10 Summarising Bundled Beta Values in Caudate
# caudate parcel ids from the incentive integration analyses
parcels.caudate<-c(355,356)
# Looping over 'caudate' parcels and create data frame of combining effects from both parcels
for (pid in parcels.caudate) { # pid<-355
# subset data for each parcel
# blocked/sustained effects data
tmp.block <- betas.block.sub %>% filter(Liquid!="Baseline", parcel.id==pid) %>%
dplyr::rename(value.block = value) %>%
mutate(Liquid = as.character(Liquid))
# bundled effects & behavioral data for cue related activity
tmp <- betas.ev.sub %>% filter(condition!="Noreward", parcel.id==pid,
(tentplot>=TR_cue[1] & tentplot<=TR_cue[2])) %>%
group_by(subID,Liquid,Reward) %>%
summarise(mean_value = mean(value)) %>% dplyr::rename(mean_value.ev = mean_value) %>%
mutate(liqCode = as.numeric(as.character(factor(Liquid,levels=c("Saltwater","Neutral","Juice"),labels=c(-1,0,1)))),
moneyCode = as.numeric(as.character(factor(Reward,levels=c("1","2","4"),labels=c(-1,0,1))))) %>%
left_join(y = tmp.block, by = c("subID","Liquid")) %>%
mutate(value.bundle = mean_value.ev + value.block) %>%
inner_join(y = incentive9.means, by = c("subID","Liquid","Reward","liqCode","moneyCode")) %>%
inner_join(y = selfreport.motivation, by = c("subID","Liquid","Reward","money")) %>%
inner_join(y = selfreport.liking, by = c("subID","Liquid","Reward","money")) %>%
inner_join(y = selfreport.perform, by = c("subID","Liquid","Reward","money"))
#tmp$liqCode <- as.numeric(levels(tmp$liqCode)[tmp$liqCode])
#tmp$moneyCode <- as.numeric(levels(tmp$moneyCode)[tmp$moneyCode])
if (pid == parcels.caudate[1]) {
data.caudate<-tmp
} else {
data.caudate<- rbind(data.caudate,tmp)
}
}`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)
`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)# Summarise the bundled values across the accumbens parcels 350 and 351
data.caudate.sum <- data.caudate %>% group_by(subID, Liquid, Reward) %>%
summarise(value.bundle.mean = mean(value.bundle), meanRR = mean(meanRR),
meanERR = mean(meanRR), meanACC = mean(meanACC),
motive_rating = mean(motive_rating), liking_rating = mean(liking_rating),
liqCode = mean(liqCode), moneyCode = mean(moneyCode)) %>%
ungroup(Liquid,Reward) %>% group_by(subID) %>%
mutate(mot9Code = liqCode + moneyCode,
meanRR.normed = scale(meanRR, center = TRUE, scale = TRUE),
meanRR.centered = scale(meanRR, center = TRUE, scale = FALSE),
motive_rating.normed = scale(motive_rating, center = TRUE, scale = TRUE),
motive_rating.centered = scale(motive_rating, center = TRUE, scale = FALSE),
value.bundle.mean.normed = scale(value.bundle.mean, center = TRUE, scale = TRUE),
value.bundle.mean.centered = scale(value.bundle.mean, center = TRUE, scale = TRUE),
value.bundle.mean.mot9 = value.bundle.mean*mot9Code,
meanRR.mot9 = meanRR*mot9Code,
motive_rating.mot9 = motive_rating*mot9Code)`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)data.caudate.sum$motive_rating.normed[is.nan(data.caudate.sum$motive_rating.normed)]<-0
# Summarise by subject
data.caudate.sum.sub <- data.caudate %>% group_by(subID) %>%
summarise(value.bundle.mean = mean(value.bundle), meanRR = mean(meanRR),
meanERR = mean(meanRR), meanACC = mean(meanACC),
motive_rating = mean(motive_rating), liking_rating = mean(liking_rating),
liqCode = mean(liqCode), moneyCode = mean(moneyCode)) %>%
#ungroup(Liquid,Reward) %>% group_by(subID) %>%
mutate(meanRR.normed = scale(meanRR, center = TRUE, scale = TRUE),
motive_rating.normed = scale(motive_rating, center = TRUE, scale = TRUE),
value.bundle.mean.normed = scale(value.bundle.mean, center = TRUE, scale = TRUE))`summarise()` ungrouping output (override with `.groups` argument)1.11 Summarising Bundled Beta Values in Putamen
# 'putamen' parcel ids from the incentive integration analyses
parcels.putamen<-c(361,362) # putamen
# Looping over 'putamen' parcels and create data frame of combining effects from both parcels
for (pid in parcels.putamen) { # pid<-361
# subset data for each parcel
# blocked/sustained effects data
tmp.block <- betas.block.sub %>% filter(Liquid!="Baseline", parcel.id==pid) %>%
dplyr::rename(value.block = value) %>%
mutate(Liquid = as.character(Liquid))
# bundled effects & behavioral data for cue related activity
tmp <- betas.ev.sub %>% filter(condition!="Noreward", parcel.id==pid,
(tentplot>=TR_cue[1] & tentplot<=TR_cue[2])) %>%
group_by(subID,Liquid,Reward) %>%
summarise(mean_value = mean(value)) %>% dplyr::rename(mean_value.ev = mean_value) %>%
mutate(liqCode = as.numeric(as.character(factor(Liquid,levels=c("Saltwater","Neutral","Juice"),labels=c(-1,0,1)))),
moneyCode = as.numeric(as.character(factor(Reward,levels=c("1","2","4"),labels=c(-1,0,1))))) %>%
left_join(y = tmp.block, by = c("subID","Liquid")) %>%
mutate(value.bundle = mean_value.ev + value.block) %>%
inner_join(y = incentive9.means, by = c("subID","Liquid","Reward","liqCode","moneyCode")) %>%
inner_join(y = selfreport.motivation, by = c("subID","Liquid","Reward","money")) %>%
inner_join(y = selfreport.liking, by = c("subID","Liquid","Reward","money")) %>%
inner_join(y = selfreport.perform, by = c("subID","Liquid","Reward","money"))
#tmp$liqCode <- as.numeric(levels(tmp$liqCode)[tmp$liqCode])
#tmp$moneyCode <- as.numeric(levels(tmp$moneyCode)[tmp$moneyCode])
if (pid == parcels.putamen[1]) {
data.putamen<-tmp
} else {
data.putamen<- rbind(data.putamen,tmp)
}
}`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)
`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)# Summarise the bundled values across the accumbens parcels 350 and 351
data.putamen.sum <- data.putamen %>% group_by(subID, Liquid, Reward) %>%
summarise(value.bundle.mean = mean(value.bundle), meanRR = mean(meanRR),
meanERR = mean(meanRR), meanACC = mean(meanACC),
motive_rating = mean(motive_rating), liking_rating = mean(liking_rating),
liqCode = mean(liqCode), moneyCode = mean(moneyCode)) %>%
ungroup(Liquid,Reward) %>% group_by(subID) %>%
mutate(mot9Code = liqCode + moneyCode,
meanRR.normed = scale(meanRR, center = TRUE, scale = TRUE),
meanRR.centered = scale(meanRR, center = TRUE, scale = FALSE),
motive_rating.normed = scale(motive_rating, center = TRUE, scale = TRUE),
motive_rating.centered = scale(motive_rating, center = TRUE, scale = FALSE),
value.bundle.mean.normed = scale(value.bundle.mean, center = TRUE, scale = TRUE),
value.bundle.mean.centered = scale(value.bundle.mean, center = TRUE, scale = TRUE),
value.bundle.mean.mot9 = value.bundle.mean*mot9Code,
meanRR.mot9 = meanRR*mot9Code,
motive_rating.mot9 = motive_rating*mot9Code)`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)data.putamen.sum$motive_rating.normed[is.nan(data.putamen.sum$motive_rating.normed)]<-0
# Summarise by subject
data.putamen.sum.sub <- data.putamen %>% group_by(subID) %>%
summarise(value.bundle.mean = mean(value.bundle), meanRR = mean(meanRR),
meanERR = mean(meanRR), meanACC = mean(meanACC),
motive_rating = mean(motive_rating), liking_rating = mean(liking_rating),
liqCode = mean(liqCode), moneyCode = mean(moneyCode)) %>%
#ungroup(Liquid,Reward) %>% group_by(subID) %>%
mutate(meanRR.normed = scale(meanRR, center = TRUE, scale = TRUE),
motive_rating.normed = scale(motive_rating, center = TRUE, scale = TRUE),
value.bundle.mean.normed = scale(value.bundle.mean, center = TRUE, scale = TRUE))`summarise()` ungrouping output (override with `.groups` argument)1.12 Summarising Bundled Beta Values in Anterior Insula
# 'anterior insula' parcel ids from the incentive integration analyses
parcels.insula<-c(99,101,143,306,340) # bilateral insula
#parcels.insula_left<-c(99,101,143) # left insula
#parcels.insula_right<-c(306,340) # right insula
#insula_hemispheres<-c("left","right")
# Looping over 'insula' parcels and create data frame of combining effects from both parcels
for (pid in parcels.insula) { # pid<-361
# subset data for each parcel
# blocked/sustained effects data
tmp.block <- betas.block.sch %>% filter(Liquid!="Baseline", parcel.id==pid) %>%
dplyr::rename(value.block = value) %>%
mutate(Liquid = as.character(Liquid))
# bundled effects & behavioral data for cue related activity
tmp <- betas.ev.sch %>% filter(condition!="Noreward", parcel.id==pid,
(tentplot>=TR_cue[1] & tentplot<=TR_cue[2])) %>%
group_by(subID,Liquid,Reward) %>%
summarise(mean_value = mean(value)) %>% dplyr::rename(mean_value.ev = mean_value) %>%
mutate(liqCode = as.numeric(as.character(factor(Liquid,levels=c("Saltwater","Neutral","Juice"),labels=c(-1,0,1)))),
moneyCode = as.numeric(as.character(factor(Reward,levels=c("1","2","4"),labels=c(-1,0,1))))) %>%
left_join(y = tmp.block, by = c("subID","Liquid")) %>%
mutate(value.bundle = mean_value.ev + value.block) %>%
inner_join(y = incentive9.means, by = c("subID","Liquid","Reward","liqCode","moneyCode")) %>%
inner_join(y = selfreport.motivation, by = c("subID","Liquid","Reward","money")) %>%
inner_join(y = selfreport.liking, by = c("subID","Liquid","Reward","money")) %>%
inner_join(y = selfreport.perform, by = c("subID","Liquid","Reward","money"))
#tmp$liqCode <- as.numeric(levels(tmp$liqCode)[tmp$liqCode])
#tmp$moneyCode <- as.numeric(levels(tmp$moneyCode)[tmp$moneyCode])
if (pid == parcels.insula[1]) {
data.insula<-tmp
} else {
data.insula<- rbind(data.insula,tmp)
}
}`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)
`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)
`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)
`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)
`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)# Summarise the bundled values across the accumbens parcels 350 and 351
data.insula.sum <- data.insula %>% group_by(subID, Liquid, Reward) %>%
summarise(value.bundle.mean = mean(value.bundle), meanRR = mean(meanRR),
meanERR = mean(meanRR), meanACC = mean(meanACC),
motive_rating = mean(motive_rating), liking_rating = mean(liking_rating),
liqCode = mean(liqCode), moneyCode = mean(moneyCode)) %>%
ungroup(Liquid,Reward) %>% group_by(subID) %>%
mutate(mot9Code = liqCode + moneyCode,
meanRR.normed = scale(meanRR, center = TRUE, scale = TRUE),
meanRR.centered = scale(meanRR, center = TRUE, scale = FALSE),
motive_rating.normed = scale(motive_rating, center = TRUE, scale = TRUE),
motive_rating.centered = scale(motive_rating, center = TRUE, scale = FALSE),
value.bundle.mean.normed = scale(value.bundle.mean, center = TRUE, scale = TRUE),
value.bundle.mean.centered = scale(value.bundle.mean, center = TRUE, scale = TRUE),
value.bundle.mean.mot9 = value.bundle.mean*mot9Code,
meanRR.mot9 = meanRR*mot9Code,
motive_rating.mot9 = motive_rating*mot9Code)`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)data.insula.sum$motive_rating.normed[is.nan(data.insula.sum$motive_rating.normed)]<-0
# Summarise by subject
data.insula.sum.sub <- data.insula %>% group_by(subID) %>%
summarise(value.bundle.mean = mean(value.bundle), meanRR = mean(meanRR),
meanERR = mean(meanRR), meanACC = mean(meanACC),
motive_rating = mean(motive_rating), liking_rating = mean(liking_rating),
liqCode = mean(liqCode), moneyCode = mean(moneyCode)) %>%
#ungroup(Liquid,Reward) %>% group_by(subID) %>%
mutate(meanRR.normed = scale(meanRR, center = TRUE, scale = TRUE),
motive_rating.normed = scale(motive_rating, center = TRUE, scale = TRUE),
value.bundle.mean.normed = scale(value.bundle.mean, center = TRUE, scale = TRUE))`summarise()` ungrouping output (override with `.groups` argument)2 Motivational Incentive Integration Effects: Reward Rate and Self-Reported Motivation Ratings
binom.test(x = 117, n = 288, p = .05, alternative = "two.sided")
Exact binomial test
data: 117 and 288
number of successes = 117, number of trials = 288, p-value < 2.2e-16
alternative hypothesis: true probability of success is not equal to 0.05
95 percent confidence interval:
0.3490173 0.4654323
sample estimates:
probability of success
0.40625 2.1 Extended Data Figure 1-2: RT between baseline and incentive block
cond.RT <- data.behavior %>% dplyr::select(subID,runID,trial,condition=runType,RT,ACC,switch_trial) %>%
filter(ACC==1) %>% mutate(RT_log = log(RT))
RT.sum=summarySEwithin2(data=cond.RT, measurevar = "RT", withinvars = c("condition"), idvar = "subID")Automatically converting the following non-factors to factors: condition
Attaching package: 'data.table'The following objects are masked from 'package:dplyr':
between, first, lastThe following object is masked from 'package:purrr':
transposepandoc.table(RT.sum)
--------------------------------------------------------------
condition N RT RTNormed sd se ci
----------- ------- ------- ---------- ------- ------- -------
baseline 3467 857.8 858.7 374.9 6.367 12.48
incentive 11443 641.7 641.4 262 2.449 4.8
--------------------------------------------------------------p.RT.0<-ggplot(data = RT.sum, aes(x=condition, y = RT, fill=condition)) +
geom_bar(position=position_dodge(width=0.8), color="black",
stat="identity", width=0.8) +
geom_errorbar(position=position_dodge(width=0.8),
aes(ymin=RT-ci, ymax=RT+ci), width=.2) +
xlab("Block") + ylab("Response Time (ms)") +
labs(fill = "Block") +
coord_cartesian(ylim=c(500,1000)) +
scale_fill_brewer(palette="Set1") +
theme_classic() +
theme(#plot.title=element_text(size=22,face="bold", vjust=2),
axis.title=element_text(size=12,face = "bold"),
axis.text=element_text(size=14),
legend.title=element_blank(),
legend.position="none",
legend.box.background = element_rect(colour = "black"),
strip.text.x = element_text(size = 12))
p.RT.0
2.2 Extended Data Paired T-test: RT between baseline and incentive block
# summarizing to have one datapoint per subject per condition.
cond.RT.sum <- cond.RT %>%group_by(subID,condition) %>%
dplyr::summarise(meanRT = mean(RT), meanRT_log = mean(RT_log))`summarise()` regrouping output by 'subID' (override with `.groups` argument)# baseline vs incentive block
t.test(x = subset(cond.RT.sum, condition=="baseline")$meanRT_log,
y = subset(cond.RT.sum, condition=="incentive")$meanRT_log,
paired = TRUE)
Paired t-test
data: subset(cond.RT.sum, condition == "baseline")$meanRT_log and subset(cond.RT.sum, condition == "incentive")$meanRT_log
t = 16.672, df = 46, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.2573245 0.3279937
sample estimates:
mean of the differences
0.2926591 # calculate cohen's d
cohen.d(subset(cond.RT.sum, condition=="baseline")$meanRT_log,
subset(cond.RT.sum, condition=="incentive")$meanRT_log,
na.rm = TRUE, pooled = TRUE, paired = TRUE, hedges.correction = FALSE)
Cohen's d
d estimate: 1.571705 (large)
95 percent confidence interval:
lower upper
1.291783 1.851628 2.3 Extended Data Figure 1-2: Accuracy between baseline and incentive block
# Compute averages for baseline and incentive blocks across all subjects
cond.ACC <- data.behavior %>% dplyr::select(subID,runID,trial,condition=runType,RT,ACC,switch_trial)
ACC.sum=summarySEwithin2(data=cond.ACC, measurevar = "ACC", withinvars = c("condition"), idvar = "subID")Automatically converting the following non-factors to factors: conditionp.ACC.0<-ggplot(data = ACC.sum, aes(x=condition, y = ACC, fill=condition)) +
geom_bar(position=position_dodge(width=0.8), color="black",
stat="identity", width = 0.8) +
geom_errorbar(position=position_dodge(width=0.8),
aes(ymin=ACC-ci, ymax=ACC+ci), width=.2) +
xlab("Block") + ylab("Accuracy") +
labs(fill = "Block") +
coord_cartesian(ylim=c(.7,1)) +
scale_fill_brewer(palette="Set1") +
theme_classic() +
theme(#plot.title=element_text(size=22,face="bold", vjust=2),
axis.title=element_text(size=12,face = "bold"),
axis.text=element_text(size=14),
legend.title=element_blank(),
legend.position="none",
legend.box.background = element_rect(colour = "black"),
strip.text.x = element_text(size = 12))
p.ACC.0
2.4 Extended Data Paired T-test: Accuracy between baseline and incentive block
# summarizing to have one datapoint per subject per condition.
cond.ACC.sum <- cond.ACC %>% group_by(subID,condition) %>%
dplyr::summarise(meanACC = mean(ACC, na.rm=TRUE))`summarise()` regrouping output by 'subID' (override with `.groups` argument)# baseline vs incentive block
t.test(x = subset(cond.ACC.sum, condition=="baseline")$meanACC,
y = subset(cond.ACC.sum, condition=="incentive")$meanACC,
paired = TRUE)
Paired t-test
data: subset(cond.ACC.sum, condition == "baseline")$meanACC and subset(cond.ACC.sum, condition == "incentive")$meanACC
t = 7.81, df = 46, p-value = 5.645e-10
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.05551917 0.09407461
sample estimates:
mean of the differences
0.07479689 # calculate cohen's d
cohen.d(subset(cond.ACC.sum, condition=="baseline")$meanACC,
subset(cond.ACC.sum, condition=="incentive")$meanACC,
na.rm = TRUE, pooled = TRUE, paired = TRUE, hedges.correction = FALSE)
Cohen's d
d estimate: 1.111597 (large)
95 percent confidence interval:
lower upper
0.7520448 1.4711495 2.5 Extended Data Figure 1-2: RT switch costs between baseline and incentive block
cond.RT.switch <- cond.RT %>%
group_by(subID,condition,taskswitch=switch_trial) %>% filter(trial!=1, ACC==1) %>%
summarise(meanRT = mean(RT)) %>% spread(key = taskswitch, value = meanRT) %>%
rename(taskswitch="switch",taskrepeat="repeat") %>%
mutate(switchcost = taskswitch-taskrepeat)`summarise()` regrouping output by 'subID', 'condition' (override with `.groups` argument)SW.sum.0<-cond.RT.switch %>% group_by(condition) %>%
summarise(SW.mean = mean(switchcost), SW.sd = sd(switchcost))`summarise()` ungrouping output (override with `.groups` argument)SW.sum<-summarySEwithin2(data=cond.RT.switch, measurevar = "switchcost",
withinvars = c("condition"), idvar = "subID")Automatically converting the following non-factors to factors: conditionpandoc.table(SW.sum)
------------------------------------------------------------------------
condition N switchcost switchcostNormed sd se ci
----------- ---- ------------ ------------------ ------- ------- -------
baseline 47 52.95 52.95 47.16 6.879 13.85
incentive 47 22.19 22.19 47.16 6.879 13.85
------------------------------------------------------------------------p.SW.1<-ggplot(data = SW.sum, aes(x=condition, y = switchcost, fill=condition)) +
geom_bar(position=position_dodge(width=0.8), color="black",
stat="identity", width=0.8) +
geom_errorbar(position=position_dodge(width=0.8),
aes(ymin=switchcost-ci, ymax=switchcost+ci), width=.2) +
xlab("Block") + ylab("RT switch cost (ms)") +
labs(fill = "Block") + scale_fill_brewer(palette="Set1") +
theme_classic() +
theme(#plot.title=element_text(size=22,face="bold", vjust=2),
axis.title=element_text(size=12,face = "bold"),
axis.text=element_text(size=14),
legend.position="none",
legend.box.background = element_rect(colour = "black"),
strip.text.x = element_text(size = 12))
p.SW.1
2.6 Extended Data Paired T-test: RT switch costs between baseline and incentive block
# baseline vs incentive block
t.test(x = subset(cond.RT.switch, condition=="baseline")$switchcost,
y = subset(cond.RT.switch, condition=="incentive")$switchcost,
paired = TRUE)
Paired t-test
data: subset(cond.RT.switch, condition == "baseline")$switchcost and subset(cond.RT.switch, condition == "incentive")$switchcost
t = 3.1615, df = 46, p-value = 0.002777
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
11.17366 50.33584
sample estimates:
mean of the differences
30.75475 # calculate cohen's d
cohen.d(subset(cond.RT.switch, condition=="baseline")$switchcost,
subset(cond.RT.switch, condition=="incentive")$switchcost,
na.rm = TRUE, pooled = TRUE, paired = TRUE, hedges.correction = FALSE)
Cohen's d
d estimate: 0.5636807 (medium)
95 percent confidence interval:
lower upper
0.1824817 0.9448797 2.7 Plot RT by switch and non-switch trials
RT.switch.sum=summarySEwithin2(data=incentive9RT.means.switch, measurevar = "meanRT",
withinvars = c("Liquid","money","switch_trial"), idvar = "subID")Automatically converting the following non-factors to factors: money, switch_trialp.RT.switch.bar<-ggplot(data = RT.switch.sum, aes(x = money, y = meanRT, fill=Liquid)) +
geom_bar(position=position_dodge(width=0.8), color="black", stat="identity", width=0.8) +
geom_errorbar(position=position_dodge(width=0.8),
aes(ymin=meanRT-se, ymax=meanRT+se), width=.2) +
facet_grid(.~switch_trial) +
xlab("Money") + ylab("RT (ms)") +
labs(fill = "Liquid") +
scale_fill_brewer(palette="Set2") +
coord_cartesian(ylim=c(600,700)) +
#scale_fill_discrete(name="Monetary Reward") +
theme_classic() +
theme(#plot.title=element_text(size=22,face="bold", vjust=2),
axis.title=element_text(size=14,face = "bold"),
axis.text=element_text(size=14),
legend.position="right",
legend.title=element_text(size=14, face = "bold"),
legend.text = element_text(size=14),
strip.text.x = element_text(size = 12))
p.RT.switch.bar
m.SW.2<-lme4::lmer(formula = meanRT ~ switchcode * liqcode * moneycode + (1|subID),
data = incentive9RT.means.switch)
tab_model(m.SW.2)| Â | meanRT | ||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 629.48 | 592.72 – 666.24 | <0.001 |
| switchcode | 22.56 | 13.67 – 31.45 | <0.001 |
| liqcode | -5.52 | -13.22 – 2.18 | 0.160 |
| moneycode | -14.31 | -22.00 – -6.61 | <0.001 |
| switchcode * liqcode | 7.07 | -3.82 – 17.96 | 0.203 |
| switchcode * moneycode | -3.36 | -14.25 – 7.52 | 0.545 |
| liqcode * moneycode | 0.54 | -8.89 – 9.97 | 0.910 |
|
(switchcode * liqcode) * moneycode |
10.23 | -3.10 – 23.57 | 0.133 |
| Random Effects | |||
| σ2 | 4258.48 | ||
| Ï„00 subID | 15709.48 | ||
| ICC | 0.79 | ||
| N subID | 46 | ||
| Observations | 828 | ||
| Marginal R2 / Conditional R2 | 0.017 / 0.790 | ||
2.8 Extended Data Save Figure 1-2: Combine RT, Accuracy, and RT switch cost plots
p.figure2<-cowplot::plot_grid(p.RT.0, p.ACC.0, p.SW.1, labels = c("a","b","c"), align = "h", ncol = 3)
ggsave(filename = "Supp_Figure2.eps", plot = p.figure2,
device="eps", dpi = 300,
path = figure.path, width = 9, height = 3, scale = 1)2.9 Reward Rate by Money and Liquid Incentives
m.RR<-lme4::lmer(formula = meanRR ~ moneyCode * liqCode + (1|subID),
data = data.dACC.sum)
m.RR.JvN<-lme4::lmer(formula = meanRR ~ moneyCode * liqCode + (1|subID),
data = data.dACC.sum, subset = Liquid!="Saltwater")
m.RR.SvN<-lme4::lmer(formula = meanRR ~ moneyCode * liqCode + (1|subID),
data = data.dACC.sum, subset = Liquid!="Juice")
tab_model(m.RR,m.RR.JvN,m.RR.SvN,
title = "Reward Rate by Experimental Factors",
dv.labels = c("Reward Rate","Juice v. Neutral","Saltwater v. Neutral"),
pred.labels = c("Intercept","Money", "Liquid","Money:Liquid"),
show.re.var = TRUE, show.stat = TRUE, show.icc = FALSE, show.ci = FALSE)| Â | Reward Rate | Juice v. Neutral | Saltwater v. Neutral | ||||||
|---|---|---|---|---|---|---|---|---|---|
| Predictors | Estimates | Statistic | p | Estimates | Statistic | p | Estimates | Statistic | p |
| Intercept | 0.69 | 43.02 | <0.001 | 0.70 | 37.76 | <0.001 | 0.70 | 40.39 | <0.001 |
| Money | 0.03 | 4.41 | <0.001 | 0.02 | 2.63 | 0.009 | 0.02 | 2.12 | 0.034 |
| Liquid | 0.01 | 2.29 | 0.022 | 0.00 | 0.08 | 0.936 | 0.03 | 2.11 | 0.034 |
| Money:Liquid | -0.00 | -0.21 | 0.830 | 0.00 | 0.27 | 0.786 | -0.01 | -0.42 | 0.673 |
| Random Effects | |||||||||
| σ2 | 0.01 | 0.01 | 0.01 | ||||||
| Ï„00 | 0.01 subID | 0.01 subID | 0.01 subID | ||||||
| N | 46 subID | 46 subID | 46 subID | ||||||
| Observations | 414 | 276 | 276 | ||||||
| Marginal R2 / Conditional R2 | 0.028 / 0.533 | 0.020 / 0.659 | 0.030 / 0.495 | ||||||
2.10 RT and Accuracy by Money and Liquid Incentives
m.RT<-lme4::lmer(formula = meanRT ~ moneyCode * liqCode + (1|subID),
data = data.dACC.sum)
m.ACC<-lme4::lmer(formula = meanACC ~ moneyCode * liqCode + (1|subID),
data = data.dACC.sum)2.11 Extended Data Table 1-2: RT and Accuracy by Money and Liquid Incentives
tab_model(m.RT, m.ACC,
dv.labels = c("Response Time","Accuracy"),
pred.labels = c("Intercept","Money", "Liquid","Money:Liquid"),
show.re.var = FALSE, show.stat = TRUE, show.icc = FALSE, show.ci = FALSE)| Â | Response Time | Accuracy | ||||
|---|---|---|---|---|---|---|
| Predictors | Estimates | Statistic | p | Estimates | Statistic | p |
| Intercept | 641.67 | 34.34 | <0.001 | 0.86 | 88.02 | <0.001 |
| Money | -16.25 | -4.93 | <0.001 | 0.01 | 1.36 | 0.175 |
| Liquid | -1.83 | -0.56 | 0.579 | 0.01 | 1.95 | 0.051 |
| Money:Liquid | 5.14 | 1.27 | 0.202 | 0.01 | 1.81 | 0.071 |
| N | 46 subID | 46 subID | ||||
| Observations | 414 | 414 | ||||
| Marginal R2 / Conditional R2 | 0.010 / 0.842 | 0.012 / 0.450 | ||||
2.12 Extended Table 1-3: Multiple Regression Between Reward Rate and Motivation Ratings
m.motive<-lme4::lmer(formula = motive_rating ~ moneyCode * liqCode + (1|subID),
data = data.dACC.sum)
m.RR.motive<-lme4::lmer(formula = meanRR ~ moneyCode * liqCode + motive_rating + (1|subID),
data = data.dACC.sum)
tab_model(m.motive,m.RR, m.RR.motive,
title = "Hierarchical Regression with Motivation Ratings",
dv.labels = c("Motivation Ratings","Reward Rate", "Reward Rate"),
pred.labels = c("Intercept","Money", "Liquid","Money:Liquid","Motivation Ratings"),
show.re.var = FALSE, show.stat = TRUE, show.icc = FALSE, show.ci = FALSE)| Â | Motivation Ratings | Reward Rate | Reward Rate | ||||||
|---|---|---|---|---|---|---|---|---|---|
| Predictors | Estimates | Statistic | p | Estimates | Statistic | p | Estimates | Statistic | p |
| Intercept | 4.93 | 42.09 | <0.001 | 0.69 | 43.02 | <0.001 | 0.61 | 23.71 | <0.001 |
| Money | 0.69 | 9.31 | <0.001 | 0.03 | 4.41 | <0.001 | 0.01 | 2.15 | 0.031 |
| Liquid | 0.75 | 10.04 | <0.001 | 0.01 | 2.29 | 0.022 | 0.00 | 0.04 | 0.967 |
| Money:Liquid | 0.02 | 0.18 | 0.858 | -0.00 | -0.21 | 0.830 | -0.00 | -0.26 | 0.795 |
| Motivation Ratings | 0.02 | 4.46 | <0.001 | ||||||
| N | 46 subID | 46 subID | 46 subID | ||||||
| Observations | 414 | 414 | 414 | ||||||
| Marginal R2 / Conditional R2 | 0.259 / 0.431 | 0.028 / 0.533 | 0.059 / 0.549 | ||||||
2.13 Hierarchical Regression with Reward Rate and Motivation Ratings
anova(m.RR,m.RR.motive)refitting model(s) with ML (instead of REML)Data: data.dACC.sum
Models:
m.RR: meanRR ~ moneyCode * liqCode + (1 | subID)
m.RR.motive: meanRR ~ moneyCode * liqCode + motive_rating + (1 | subID)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
m.RR 6 -612.69 -588.53 312.34 -624.69
m.RR.motive 7 -630.38 -602.20 322.19 -644.38 19.698 1 9.072e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 12.14 Reward Rate by Motivation Ratings By Subject (Scatterplot)
# Reward Rate by Self Report Motivation Ratings by Subject
# RColorBrewer::brewer.pal(8, "Set1") # to get the hex colors for set 1
cor.RR.motive<-cor.test(x = data.dACC.sum.sub$motive_rating.normed, y = data.dACC.sum.sub$meanRR.normed)
Rsquared=round(cor.RR.motive$estimate^2,3)
textlab1 <- paste("R^2 == ", Rsquared)
textlab2 <- paste("t(44)= ",round(cor.RR.motive$statistic,3),", p = ",round(cor.RR.motive$p.value,2), sep="")
p.linecolor<-brewer.pal(n=8, name = 'Set1')
p.RR.motive.bysub<-ggplot(data = data.dACC.sum.sub, aes(x = motive_rating.normed, y = meanRR.normed)) +
geom_point(size = 2) +
geom_smooth(method="lm", se = TRUE, col = p.linecolor[2]) +
annotate("text", x = .15, y = -2.0, label = textlab1, color=p.linecolor[2], size = 4, parse=TRUE, hjust=0) +
annotate("text", x = .15, y = -2.3, label = textlab2, color=p.linecolor[2], size = 4, parse=FALSE, hjust=0) +
xlab("Self-Report Motivation Ratings") + ylab("Reward Rate") +
#ggtitle("Average Values by Subject") +
scale_color_brewer(palette="Set1") +
theme_classic() +
theme(panel.grid.minor=element_blank(),
plot.background = element_rect(fill = "transparent",colour = NA),
rect = element_rect(fill = "transparent"),
axis.title=element_text(size=12,face = "bold"),
axis.text=element_text(size=16),
legend.position="none",
strip.text.x = element_text(size = 16))
p.RR.motive.bysub`geom_smooth()` using formula 'y ~ x'
2.15 Reward Rate by Motivation Ratings by Monetary Reward and Liquid (Scatterplot)
#p.linecolor<-brewer.pal(n=9, name = "Spectral")
p.RR.motive<-ggplot(data = data.dACC.sum,
aes(x = motive_rating.normed, y = meanRR.normed, color=Liquid)) +
geom_point(size = 2) +
stat_smooth(method="glm", se = FALSE) +
facet_grid(.~MoneyLabel) +
#annotate("text", x = .35, y = -2.0, label = textlab1, color=p.linecolor[4], size = 5, parse=FALSE, hjust=0) +
#annotate("text", x = .35, y = -2.3, label = textlab2, color=p.linecolor[4], size = 5, parse=FALSE, hjust=0) +
labs(x = "Self Report Motivation Ratings", y = "Reward Rate", color = "Liquid") +
scale_color_brewer(palette="Set2") +
#scale_color_manual(labels = c("$", "$$","$$$$")) +
theme_bw() +
theme(#plot.title=element_text(size=22,face="bold", vjust=2),
axis.title=element_text(size=14,face = "bold"),
axis.text=element_text(size=14),
legend.position="right",
legend.title = element_text(size = 14, face = "bold"),
legend.text = element_text(size=14),
strip.text.x = element_text(size = 14))
p.RR.motive`geom_smooth()` using formula 'y ~ x'
2.16 Motivation Ratings by Monetary Reward and Liquid (Barplot)
RR.sum=summarySEwithin2(data=data.dACC.sum, measurevar = "motive_rating",
withinvars = c("Liquid","MoneyLabel"), idvar = "subID")Automatically converting the following non-factors to factors: Liquidp.RR.motive.bar<-ggplot(data = RR.sum, aes(x = MoneyLabel, y = motive_rating, fill=Liquid)) +
geom_bar(position=position_dodge(width=0.8), color="black", stat="identity", width=0.8) +
geom_errorbar(position=position_dodge(width=0.8),
aes(ymin=motive_rating-se, ymax=motive_rating+se), width=.2) +
xlab("Money") + ylab("Self-Report Motivation Ratings") +
labs(fill = "Liquid") +
#geom_hline(yintercept = .4, linetype="dashed") +
scale_fill_brewer(palette="Set2") +
coord_cartesian(ylim=c(1,7)) +
#scale_fill_discrete(name="Monetary Reward") +
theme_classic() +
theme(#plot.title=element_text(size=22,face="bold", vjust=2),
axis.title=element_text(size=14,face = "bold"),
axis.text=element_text(size=14),
legend.position="right",
legend.title=element_text(size=14, face = "bold"),
legend.text = element_text(size=14),
strip.text.x = element_text(size = 12))
p.RR.motive.bar
2.17 Reward Rate by Monetary Reward and Liquid (Barplot)
RR.sum=summarySEwithin2(data=data.dACC.sum, measurevar = "meanRR",
withinvars = c("Liquid","MoneyLabel"), idvar = "subID")Automatically converting the following non-factors to factors: Liquidp.RR.mot9<-ggplot(RR.sum, aes(x=MoneyLabel, y=meanRR, fill=Liquid)) +
geom_bar(position=position_dodge(width=0.8), color="black",
stat="identity", width=0.8) +
geom_errorbar(position=position_dodge(width=0.8),
aes(ymin=meanRR-se, ymax=meanRR+se), width=.2) +
xlab("Money") + ylab("Reward Rate") +
labs(fill = "Liquid") +
geom_hline(yintercept = .4, linetype="dashed") +
scale_fill_brewer(palette="Set2") +
coord_cartesian(ylim=c(.5,.75)) +
#scale_fill_discrete(name="Monetary Reward") +
theme_classic() +
theme(#plot.title=element_text(size=22,face="bold", vjust=2),
axis.title=element_text(size=14,face = "bold"),
axis.text=element_text(size=14),
legend.position="none",
legend.title=element_text(size=14, face = "bold"),
legend.text = element_text(size=14),
strip.text.x = element_text(size = 12))
p.RR.mot9
2.18 Manuscript Save Figure 1
p.figure1<-cowplot::plot_grid(p.RR.mot9, p.RR.motive.bar, p.RR.motive.bysub,
label_size = 16, ncol = 3, rel_widths = c(1,1.2,1))`geom_smooth()` using formula 'y ~ x'p.figure1
#labels = c("B","C")
ggsave(filename = "Figure1.eps", plot = p.figure1,
device="eps", dpi = 300,
path = figure.path, width = 12, height = 3, scale = 1.2)Warning in grid.Call.graphics(C_polygon, x$x, x$y, index): semi-transparency is
not supported on this device: reported only once per page2.19 Extended Data Figure 1-4: RT and Accuracy by Monetary Reward and Liquid (Barplot)
RT.sum=summarySEwithin2(data=data.dACC.sum, measurevar = "meanRT",
withinvars = c("Liquid","MoneyLabel"), idvar = "subID")Automatically converting the following non-factors to factors: Liquidp.RT.mot9<-ggplot(RT.sum, aes(x=MoneyLabel, y=meanRT, fill=Liquid)) +
geom_bar(position=position_dodge(width=0.8), color="black",
stat="identity", width=0.8) +
geom_errorbar(position=position_dodge(width=0.8),
aes(ymin=meanRT-se, ymax=meanRT+se), width=.2) +
xlab("Money") + ylab("Response Time (ms)") +
labs(fill = "Liquid") +
geom_hline(yintercept = .4, linetype="dashed") +
scale_fill_brewer(palette="Set2") +
coord_cartesian(ylim=c(600,700)) +
#scale_fill_discrete(name="Monetary Reward") +
theme_classic() +
theme(#plot.title=element_text(size=22,face="bold", vjust=2),
axis.title=element_text(size=14,face = "bold"),
axis.text=element_text(size=14),
legend.position="none",
legend.title=element_text(size=14, face = "bold"),
legend.text = element_text(size=14),
strip.text.x = element_text(size = 12))
p.RT.mot9
ACC.sum=summarySEwithin2(data=data.dACC.sum, measurevar = "meanACC",
withinvars = c("Liquid","MoneyLabel"), idvar = "subID")Automatically converting the following non-factors to factors: Liquidp.ACC.mot9<-ggplot(ACC.sum, aes(x=MoneyLabel, y=meanACC, fill=Liquid)) +
geom_bar(position=position_dodge(width=0.8), color="black",
stat="identity", width=0.8) +
geom_errorbar(position=position_dodge(width=0.8),
aes(ymin=meanACC-se, ymax=meanACC+se), width=.2) +
xlab("Money") + ylab("Accuracy") +
labs(fill = "Liquid") +
geom_hline(yintercept = .4, linetype="dashed") +
scale_fill_brewer(palette="Set2") +
coord_cartesian(ylim=c(.75,1)) +
#scale_fill_discrete(name="Monetary Reward") +
theme_classic() +
theme(#plot.title=element_text(size=22,face="bold", vjust=2),
axis.title=element_text(size=14,face = "bold"),
axis.text=element_text(size=14),
legend.position="right",
legend.title=element_text(size=14, face = "bold"),
legend.text = element_text(size=14),
strip.text.x = element_text(size = 12))
p.ACC.mot9
2.20 Extended Data Save Figure 1-4
p.figure4_sup<-cowplot::plot_grid(p.RT.mot9, p.ACC.mot9,
label_size = 16, ncol = 2, rel_widths = c(1,1.3))
p.figure4_sup
ggsave(filename = "Supplemental_Figure1-4.eps", plot = p.figure4_sup,
device="eps", dpi = 300,
path = figure.path, width = 7, height = 3, scale = 1.2)3 Dorsal Anterior Cingulate Cortex (dACC) Encodes Both Monetary and Liquid Incentives in Bundled Beta Estimates
3.1 Block Beta Values by Liquid in dACC
m.dACC.RR.block<-lme4::lmer(formula = value.block ~ moneyCode * liqCode + (1|subID),
data = data.dACC.sum)3.3 Bundled Beta Values by Money in dACC
m.dACC.RR.bundle<-lme4::lmer(formula = value.bundle.mean ~ moneyCode * liqCode + (1|subID),
data = data.dACC.sum)3.5 Manuscript Figure 2: Plotting dACC Blocked (Liquid) Effects & Cue (Money) Effects
dACC.block.sum=summarySEwithin2(data = data.dACC.sum, measurevar = "value.block",
withinvars = "Liquid", idvar = "subID")Automatically converting the following non-factors to factors: Liquidp.figure2.block<-ggplot(data = dACC.block.sum, aes(x = Liquid, y = value.block, fill = Liquid)) +
geom_bar(stat = "identity", color="black") +
geom_errorbar(position=position_dodge(width=0.8),
aes(ymin=value.block-se, ymax=value.block+se), width=.2) +
xlab("Liquid") + ylab("Block Beta Estimates") +
labs(fill = "Liquid") +
#geom_hline(yintercept = .4, linetype="dashed") +
scale_fill_brewer(palette="Set2") +
#coord_cartesian(ylim=c(.35,.8)) +
#scale_fill_discrete(name="Monetary Reward") +
theme_classic()+
theme(#plot.title=element_text(size=22,face="bold", vjust=2),
axis.title=element_text(size=14,face = "bold"),
axis.text=element_text(size=14),
legend.position="none",
legend.title=element_text(size=14, face = "bold"),
strip.text.x = element_text(size = 12))
p.figure2.block
dACC.ev.sum=summarySEwithin2(data = data.dACC.sum, measurevar = "value.ev",
withinvars = c("Liquid","MoneyLabel"), idvar = "subID")Automatically converting the following non-factors to factors: Liquidp.figure2.ev<-ggplot(data = dACC.ev.sum, aes(x = MoneyLabel, y = value.ev, fill = Liquid)) +
geom_bar(position=position_dodge(width=0.8), color="black",
stat="identity", width=0.8) +
geom_errorbar(position=position_dodge(width=0.8),
aes(ymin=value.ev-se, ymax=value.ev+se), width=.2) +
xlab("Money") + ylab("Cue Beta Estimates") +
labs(fill = "Liquid") +
#geom_hline(yintercept = .4, linetype="dashed") +
scale_fill_brewer(palette="Set2") +
#coord_cartesian(ylim=c(.35,.8)) +
#scale_fill_discrete(name="Monetary Reward") +
theme_classic() +
theme(#plot.title=element_text(size=22,face="bold", vjust=2),
axis.title=element_text(size=14,face = "bold"),
axis.text=element_text(size=14),
legend.position="right",
legend.title=element_text(size=14, face = "bold"),
legend.text = element_text(size=14),
strip.text.x = element_text(size = 12))
p.figure2.ev
3.6 Manuscript Figure 2: Plotting dACC Bundled Effects
dACC.bundle.sum=summarySEwithin2(data = data.dACC.sum, measurevar = "value.bundle.mean",
withinvars = c("Liquid","MoneyLabel"), idvar = "subID")Automatically converting the following non-factors to factors: Liquidp.figure2.bundle<-ggplot(data = dACC.bundle.sum, aes(x = MoneyLabel, y = value.bundle.mean, fill = Liquid)) +
geom_bar(position=position_dodge(width=0.8), color="black",
stat="identity", width=0.8) +
geom_errorbar(position=position_dodge(width=0.8),
aes(ymin=value.bundle.mean-se, ymax=value.bundle.mean+se), width=.2) +
xlab("Money") + ylab("Bundled Beta Estimates") +
labs(fill = "Liquid") +
#geom_hline(yintercept = .4, linetype="dashed") +
scale_fill_brewer(palette="Set2") +
#coord_cartesian(ylim=c(.35,.8)) +
#scale_fill_discrete(name="Monetary Reward") +
theme_classic() +
theme(#plot.title=element_text(size=22,face="bold", vjust=2),
axis.title=element_text(size=14,face = "bold"),
axis.text=element_text(size=14),
legend.position="right",
legend.title=element_text(size=14, face = "bold"),
legend.text = element_text(size=14),
strip.text.x = element_text(size = 12))
p.figure2.bundle
3.7 Manuscript Figure 2: Draw Tentzero function of trialwise beta estimates
#'dACC' parcel ids from the incentive integration analyses
parcels.dACC<-c(107,110,108,311,312)
parcels.vmpfc<-c(168, 169, 173, 174, 379, 380, 381, 382)
# colors.TENT<-c("firebrick1","firebrick3","firebrick4",
# "seagreen1","seagreen3","seagreen4",
# "steelblue1","steelblue3","steelblue4")
colors.TENT<-c("#62c1a3","#3A9276","#1D493B",
"#f98c61","#EC4909","#9D3106",
"#8B9FCB","#4C69A9","#2C3D63")
data.tentplot<-betas.ev.sch %>% filter(condition!="Noreward", parcel.id %in% parcels.dACC) %>%
mutate(time=tentplot*2) %>%
group_by(time,Liquid,Reward,condition) %>%
#group_by(time,Reward) %>%
dplyr::summarise(mean_value = mean(value), .groups = "keep")
# add zero padding to beginning and end of tent functions
data.tentplot2<-data.tentplot %>% #dplyr::filter(time<=4) %>%
ungroup() %>%
add_row(time = c(rep(0,9),rep(24,9)), Liquid = data.tentplot$Liquid[1:18],
Reward = data.tentplot$Reward[1:18], condition = data.tentplot$condition[1:18],
mean_value=0) %>%
mutate(MoneyLabel = factor(x = Reward, levels = c(1,2,4), labels = c("Low", "Med", "High")))
# data.tentplot<-data.tentplot %>% #dplyr::filter(time<=4) %>%
# ungroup(time) %>%
# mutate(time=c(rep(0,9),rep(24,9)), mean_value=0) %>%
# bind_rows(data.tentplot) %>%
# mutate(MoneyLabel = factor(x = Reward, levels = c(1,2,4), labels = c("Low", "Med", "High")))
p.trialbetas<-ggplot(data = data.tentplot2, mapping = aes(x = time, y = mean_value)) +
geom_rect(mapping = aes(xmin=3.5,xmax=6.5,ymin=-.2,ymax=0.4), fill = "#e0e1dd", color = "#e0e1dd") +
geom_line(aes(group=Reward, linetype = MoneyLabel, color = condition), size= .9) +
#geom_point(aes(group=Reward, color = condition)) +
facet_grid(.~Liquid) +
scale_color_manual(values=colors.TENT, guide = FALSE) +
#geom_text(x=5, y=.4, label="Cue", size = 6) +
labs(x = "Time from Cue Onset (ms)", y = "Trial Beta Estimates", linetype = "Money") +
scale_x_continuous(breaks=c(seq(0,24,4)), labels=c(seq(0,24,4)), limits=c(0,24)) +
scale_y_continuous(breaks=c(seq(-.2,.5,.2)), labels=c(seq(-.2,.5,.2)), limits=c(-.2,.5)) +
theme_classic() +
theme(panel.grid.minor=element_blank(),
plot.background = element_rect(fill = "transparent",colour = NA),
rect = element_rect(fill = "transparent"),
#plot.title=element_text(size=22,face="bold", vjust=2),
axis.title=element_text(size=15,face = "bold"),
axis.text=element_text(size=16),
strip.text.x = element_text(size = 16, face="bold"),
legend.title=element_text(size=16,face="bold"),
legend.text=element_text(size=12),
legend.position="top")
p.trialbetas
3.8 Manuscript Save Figure 2
ggsave(filename = "Figure2_block.eps", plot = p.figure2.block,
device="eps", dpi = 300,
path = figure.path, width = 5, height = 3, scale = 1.1)
ggsave(filename = "Figure2_ev.eps", plot = p.figure2.ev,
device="eps", dpi = 300,
path = figure.path, width = 7, height = 4, scale = 1.1)
ggsave(filename = "Figure2_bundle.eps", plot = p.figure2.bundle,
device="eps", dpi = 300,
path = figure.path, width = 7, height = 4, scale = 1.1)
ggsave(filename = "Figure2_trialbetas.eps", plot = p.trialbetas,
device="eps", dpi = 300,
path = figure.path, width = 10, height = 4, scale = 1.1,
bg = "transparent")4 Reward Rate and Self-Report Motivation Ratings Predicted by Dorsal ACC Bundled Betas
4.1 Reward Rate predicted by Bundled Betas
m.dACC.RR<-lme4::lmer(formula = meanRR ~ mot9Code + value.bundle.mean + (1+mot9Code|subID),
data = data.dACC.sum)
m.dACC.RR.null<-update(m.dACC.RR, formula = ~ . -value.bundle.mean )
BF_dACC_RR = BIC(m.dACC.RR.null) - BIC(m.dACC.RR)4.2 Self-Report Motivation Ratings Predicted by Dorsal ACC Bundled Betas
m.dACC.motive<-lme4::lmer(formula = motive_rating ~ mot9Code + value.bundle.mean + (1+mot9Code|subID),
data = data.dACC.sum)4.3 Plotting the Models
tab_model(m.dACC.RR, m.dACC.motive,
title = "Reward Rate by Bundled dACC Betas",
dv.labels = c("Reward Rate","Motivation Ratings"),
pred.labels = c("Intercept","Incentive Conditions","dACC Bundled Betas"),
show.re.var = FALSE, show.stat = TRUE, show.icc = FALSE)| Â | Reward Rate | Motivation Ratings | ||||||
|---|---|---|---|---|---|---|---|---|
| Predictors | Estimates | CI | Statistic | p | Estimates | CI | Statistic | p |
| Intercept | 0.68 | 0.65 – 0.71 | 40.93 | <0.001 | 4.70 | 4.45 – 4.94 | 37.43 | <0.001 |
| Incentive Conditions | 0.02 | 0.01 – 0.03 | 2.91 | 0.004 | 0.67 | 0.54 – 0.81 | 9.78 | <0.001 |
| dACC Bundled Betas | 0.09 | 0.02 – 0.17 | 2.41 | 0.016 | 1.64 | 0.83 – 2.45 | 3.96 | <0.001 |
| N | 46 subID | 46 subID | ||||||
| Observations | 414 | 414 | ||||||
| Marginal R2 / Conditional R2 | 0.040 / 0.591 | 0.294 / 0.501 | ||||||
5 Dorsal ACC Effects on Reward Rate Mediated by Self-Report Motivation Ratings (Within-Subjects Mediation Analyses)
x = dACC /newline y = Reward Rate /newline m = Motivation Ratings newline
5.1 Format the Data for Within-Subjects Mediation Analysis
data.dACC.mediation <- data.dACC.sum %>%
dplyr::select(subID,meanRR,motive_rating,value.bundle.mean,mot9Code,liqCode,moneyCode)5.2 Step 0: Motivation Conditions Predict Reward Rate (Y), dACC Betas, and Motivation Ratings
m.step0.RR<-lme4::lmer(formula = meanRR ~ mot9Code + (1+mot9Code|subID), data = data.dACC.sum)
m.step0.beta<-lme4::lmer(formula = value.bundle.mean ~ mot9Code + (1+mot9Code|subID), data = data.dACC.sum)
m.step0.motive<-lme4::lmer(formula = motive_rating ~ mot9Code + (1+mot9Code|subID), data = data.dACC.sum)
tab_model(m.step0.RR, m.step0.beta, m.step0.motive,
dv.labels = c("Reward Rate","dACC Betas", "Motivation Ratings"),
pred.labels = c("Intercept","Incentive Conditions"),
show.re.var = FALSE, show.stat = TRUE, show.icc = FALSE)| Â | Reward Rate | dACC Betas | Motivation Ratings | |||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Predictors | Estimates | CI | Statistic | p | Estimates | CI | Statistic | p | Estimates | CI | Statistic | p |
| Intercept | 0.69 | 0.66 – 0.73 | 43.02 | <0.001 | 0.14 | 0.10 – 0.18 | 6.57 | <0.001 | 4.93 | 4.70 – 5.16 | 42.09 | <0.001 |
| Incentive Conditions | 0.02 | 0.01 – 0.03 | 3.36 | 0.001 | 0.03 | 0.01 – 0.04 | 4.04 | <0.001 | 0.72 | 0.58 – 0.86 | 10.35 | <0.001 |
| N | 46 subID | 46 subID | 46 subID | |||||||||
| Observations | 414 | 414 | 414 | |||||||||
| Marginal R2 / Conditional R2 | 0.025 / 0.597 | 0.031 / 0.613 | 0.259 / 0.492 | |||||||||
5.3 Step 1: Reward Rate (Y) by dACC Betas (X)
m.dACC.step1<-lme4::lmer(formula = meanRR ~ value.bundle.mean + (1|subID),
data = data.dACC.sum)
tab_model(m.dACC.step1)| Â | meanRR | ||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 0.67 | 0.64 – 0.71 | <0.001 |
| value.bundle.mean | 0.14 | 0.07 – 0.22 | <0.001 |
| Random Effects | |||
| σ2 | 0.01 | ||
| Ï„00 subID | 0.01 | ||
| ICC | 0.49 | ||
| N subID | 46 | ||
| Observations | 414 | ||
| Marginal R2 / Conditional R2 | 0.035 / 0.510 | ||
5.4 Step 2: Motivation Ratings (M) Predicted by dACC Betas
m.dACC.step2<-lme4::lmer(formula = motive_rating ~ value.bundle.mean + (1|subID),
data = data.dACC.sum)5.5 Step 3: dACC Betas Predict Reward Rate Controling for Motivation Ratings
m.dACC.step3<-lme4::lmer(formula = meanRR ~ value.bundle.mean + motive_rating + (1|subID),
data = data.dACC.sum)5.6 Plot all 3 Models
tab_model(m.dACC.step1,m.dACC.step2,m.dACC.step3,
title = "dACC Within-Subjects Mediation Models",
dv.labels = c("Reward Rate","Motivation Ratings","Reward Rate"),
#pred.labels = c("Intercept","Money", "Liquid","Motivation Ratings","dACC Bundled Betas"),
show.re.var = FALSE, show.stat = TRUE, show.icc = FALSE, show.ci = FALSE)| Â | Reward Rate | Motivation Ratings | Reward Rate | ||||||
|---|---|---|---|---|---|---|---|---|---|
| Predictors | Estimates | Statistic | p | Estimates | Statistic | p | Estimates | Statistic | p |
| (Intercept) | 0.67 | 40.87 | <0.001 | 4.53 | 34.68 | <0.001 | 0.59 | 26.08 | <0.001 |
| value.bundle.mean | 0.14 | 3.77 | <0.001 | 2.77 | 5.98 | <0.001 | 0.08 | 2.15 | 0.032 |
| motive_rating | 0.02 | 5.53 | <0.001 | ||||||
| N | 46 subID | 46 subID | 46 subID | ||||||
| Observations | 414 | 414 | 414 | ||||||
| Marginal R2 / Conditional R2 | 0.035 / 0.510 | 0.102 / 0.231 | 0.071 / 0.545 | ||||||
5.7 Step 4: Run the mediation using bmlm package
data.dACC.mediation<-data.frame(data.dACC.mediation)
data.dACC.mediation<-isolate(d = data.dACC.mediation,by = "subID",
value = c("meanRR","motive_rating","value.bundle.mean","mot9Code"))
# Running the multilevel mediation on mean-centered data
set.seed(200)
dACC.mediation.results = mlm(d = data.dACC.mediation, id = "subID",
x = "value.bundle.mean_cw", m = "motive_rating_cw", y = "meanRR_cw",
iter = 10000)Estimating model, please wait.
SAMPLING FOR MODEL 'bmlm' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.000211 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 2.11 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 10000 [ 0%] (Warmup)
Chain 1: Iteration: 1000 / 10000 [ 10%] (Warmup)
Chain 1: Iteration: 2000 / 10000 [ 20%] (Warmup)
Chain 1: Iteration: 3000 / 10000 [ 30%] (Warmup)
Chain 1: Iteration: 4000 / 10000 [ 40%] (Warmup)
Chain 1: Iteration: 5000 / 10000 [ 50%] (Warmup)
Chain 1: Iteration: 5001 / 10000 [ 50%] (Sampling)
Chain 1: Iteration: 6000 / 10000 [ 60%] (Sampling)
Chain 1: Iteration: 7000 / 10000 [ 70%] (Sampling)
Chain 1: Iteration: 8000 / 10000 [ 80%] (Sampling)
Chain 1: Iteration: 9000 / 10000 [ 90%] (Sampling)
Chain 1: Iteration: 10000 / 10000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 16.9229 seconds (Warm-up)
Chain 1: 13.053 seconds (Sampling)
Chain 1: 29.9759 seconds (Total)
Chain 1:
SAMPLING FOR MODEL 'bmlm' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 8.9e-05 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.89 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 10000 [ 0%] (Warmup)
Chain 2: Iteration: 1000 / 10000 [ 10%] (Warmup)
Chain 2: Iteration: 2000 / 10000 [ 20%] (Warmup)
Chain 2: Iteration: 3000 / 10000 [ 30%] (Warmup)
Chain 2: Iteration: 4000 / 10000 [ 40%] (Warmup)
Chain 2: Iteration: 5000 / 10000 [ 50%] (Warmup)
Chain 2: Iteration: 5001 / 10000 [ 50%] (Sampling)
Chain 2: Iteration: 6000 / 10000 [ 60%] (Sampling)
Chain 2: Iteration: 7000 / 10000 [ 70%] (Sampling)
Chain 2: Iteration: 8000 / 10000 [ 80%] (Sampling)
Chain 2: Iteration: 9000 / 10000 [ 90%] (Sampling)
Chain 2: Iteration: 10000 / 10000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 15.6859 seconds (Warm-up)
Chain 2: 13.3309 seconds (Sampling)
Chain 2: 29.0168 seconds (Total)
Chain 2:
SAMPLING FOR MODEL 'bmlm' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 0.000106 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 1.06 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 10000 [ 0%] (Warmup)
Chain 3: Iteration: 1000 / 10000 [ 10%] (Warmup)
Chain 3: Iteration: 2000 / 10000 [ 20%] (Warmup)
Chain 3: Iteration: 3000 / 10000 [ 30%] (Warmup)
Chain 3: Iteration: 4000 / 10000 [ 40%] (Warmup)
Chain 3: Iteration: 5000 / 10000 [ 50%] (Warmup)
Chain 3: Iteration: 5001 / 10000 [ 50%] (Sampling)
Chain 3: Iteration: 6000 / 10000 [ 60%] (Sampling)
Chain 3: Iteration: 7000 / 10000 [ 70%] (Sampling)
Chain 3: Iteration: 8000 / 10000 [ 80%] (Sampling)
Chain 3: Iteration: 9000 / 10000 [ 90%] (Sampling)
Chain 3: Iteration: 10000 / 10000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 15.7006 seconds (Warm-up)
Chain 3: 13.9962 seconds (Sampling)
Chain 3: 29.6968 seconds (Total)
Chain 3:
SAMPLING FOR MODEL 'bmlm' NOW (CHAIN 4).
Chain 4:
Chain 4: Gradient evaluation took 0.000113 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 1.13 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 10000 [ 0%] (Warmup)
Chain 4: Iteration: 1000 / 10000 [ 10%] (Warmup)
Chain 4: Iteration: 2000 / 10000 [ 20%] (Warmup)
Chain 4: Iteration: 3000 / 10000 [ 30%] (Warmup)
Chain 4: Iteration: 4000 / 10000 [ 40%] (Warmup)
Chain 4: Iteration: 5000 / 10000 [ 50%] (Warmup)
Chain 4: Iteration: 5001 / 10000 [ 50%] (Sampling)
Chain 4: Iteration: 6000 / 10000 [ 60%] (Sampling)
Chain 4: Iteration: 7000 / 10000 [ 70%] (Sampling)
Chain 4: Iteration: 8000 / 10000 [ 80%] (Sampling)
Chain 4: Iteration: 9000 / 10000 [ 90%] (Sampling)
Chain 4: Iteration: 10000 / 10000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 17.3432 seconds (Warm-up)
Chain 4: 13.5487 seconds (Sampling)
Chain 4: 30.8919 seconds (Total)
Chain 4: # visualizing the path
mlm_path_plot(dACC.mediation.results, text = T)
# summarising within-subject effects
mlm.sum<-mlm_summary(dACC.mediation.results, level = .9)
mlm_pars_plot(dACC.mediation.results, type = "coef", pars = c( "cp", "c", "me","pme"))No id variables; using all as measure variablesWarning: `group_by_()` is deprecated as of dplyr 0.7.0.
Please use `group_by()` instead.
See vignette('programming') for more help
This warning is displayed once every 8 hours.
Call `lifecycle::last_warnings()` to see where this warning was generated.Warning: `summarise_()` is deprecated as of dplyr 0.7.0.
Please use `summarise()` instead.
This warning is displayed once every 8 hours.
Call `lifecycle::last_warnings()` to see where this warning was generated.
mlm_pars_plot(dACC.mediation.results, type = "hist", pars = c( "cp", "c", "me"))No id variables; using all as measure variables`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# check for convergence of mcmc traces, requires bayesplot package
mcmc_trace(as.data.frame(dACC.mediation.results), pars = c("a","b","cp","c"))
mcmc_rank_hist(as.data.frame(dACC.mediation.results), pars = c("a","b","cp","c"))
# Create a pretty table of the parameters estimates
formattable(mlm.sum)| Parameter | Mean | SE | Median | 5% | 95% | n_eff | Rhat |
|---|---|---|---|---|---|---|---|
| a | 3.23 | 0.99 | 3.23 | 1.60 | 4.84 | 16780 | 1 |
| b | 0.02 | 0.00 | 0.02 | 0.01 | 0.03 | 16523 | 1 |
| cp | 0.07 | 0.05 | 0.07 | -0.01 | 0.15 | 16338 | 1 |
| me | 0.05 | 0.03 | 0.05 | 0.01 | 0.10 | 8200 | 1 |
| c | 0.13 | 0.05 | 0.13 | 0.04 | 0.21 | 20909 | 1 |
| pme | 0.37 | 20.48 | 0.42 | 0.05 | 1.10 | 20428 | 1 |
5.8 Inverse mediation, with RR mediating effects of dACC on motivation ratings.
set.seed(500)
dACC.mediation.inverse.results = mlm(d = data.dACC.mediation, id = "subID",
x = "value.bundle.mean_cw", m = "meanRR_cw", y = "motive_rating_cw",
iter = 10000)Estimating model, please wait.
SAMPLING FOR MODEL 'bmlm' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.00011 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.1 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 10000 [ 0%] (Warmup)
Chain 1: Iteration: 1000 / 10000 [ 10%] (Warmup)
Chain 1: Iteration: 2000 / 10000 [ 20%] (Warmup)
Chain 1: Iteration: 3000 / 10000 [ 30%] (Warmup)
Chain 1: Iteration: 4000 / 10000 [ 40%] (Warmup)
Chain 1: Iteration: 5000 / 10000 [ 50%] (Warmup)
Chain 1: Iteration: 5001 / 10000 [ 50%] (Sampling)
Chain 1: Iteration: 6000 / 10000 [ 60%] (Sampling)
Chain 1: Iteration: 7000 / 10000 [ 70%] (Sampling)
Chain 1: Iteration: 8000 / 10000 [ 80%] (Sampling)
Chain 1: Iteration: 9000 / 10000 [ 90%] (Sampling)
Chain 1: Iteration: 10000 / 10000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 12.9923 seconds (Warm-up)
Chain 1: 6.2606 seconds (Sampling)
Chain 1: 19.2529 seconds (Total)
Chain 1:
SAMPLING FOR MODEL 'bmlm' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 8.8e-05 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.88 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 10000 [ 0%] (Warmup)
Chain 2: Iteration: 1000 / 10000 [ 10%] (Warmup)
Chain 2: Iteration: 2000 / 10000 [ 20%] (Warmup)
Chain 2: Iteration: 3000 / 10000 [ 30%] (Warmup)
Chain 2: Iteration: 4000 / 10000 [ 40%] (Warmup)
Chain 2: Iteration: 5000 / 10000 [ 50%] (Warmup)
Chain 2: Iteration: 5001 / 10000 [ 50%] (Sampling)
Chain 2: Iteration: 6000 / 10000 [ 60%] (Sampling)
Chain 2: Iteration: 7000 / 10000 [ 70%] (Sampling)
Chain 2: Iteration: 8000 / 10000 [ 80%] (Sampling)
Chain 2: Iteration: 9000 / 10000 [ 90%] (Sampling)
Chain 2: Iteration: 10000 / 10000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 11.9064 seconds (Warm-up)
Chain 2: 9.0943 seconds (Sampling)
Chain 2: 21.0007 seconds (Total)
Chain 2:
SAMPLING FOR MODEL 'bmlm' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 9.8e-05 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.98 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 10000 [ 0%] (Warmup)
Chain 3: Iteration: 1000 / 10000 [ 10%] (Warmup)
Chain 3: Iteration: 2000 / 10000 [ 20%] (Warmup)
Chain 3: Iteration: 3000 / 10000 [ 30%] (Warmup)
Chain 3: Iteration: 4000 / 10000 [ 40%] (Warmup)
Chain 3: Iteration: 5000 / 10000 [ 50%] (Warmup)
Chain 3: Iteration: 5001 / 10000 [ 50%] (Sampling)
Chain 3: Iteration: 6000 / 10000 [ 60%] (Sampling)
Chain 3: Iteration: 7000 / 10000 [ 70%] (Sampling)
Chain 3: Iteration: 8000 / 10000 [ 80%] (Sampling)
Chain 3: Iteration: 9000 / 10000 [ 90%] (Sampling)
Chain 3: Iteration: 10000 / 10000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 12.1626 seconds (Warm-up)
Chain 3: 12.2249 seconds (Sampling)
Chain 3: 24.3875 seconds (Total)
Chain 3:
SAMPLING FOR MODEL 'bmlm' NOW (CHAIN 4).
Chain 4:
Chain 4: Gradient evaluation took 0.000137 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 1.37 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 10000 [ 0%] (Warmup)
Chain 4: Iteration: 1000 / 10000 [ 10%] (Warmup)
Chain 4: Iteration: 2000 / 10000 [ 20%] (Warmup)
Chain 4: Iteration: 3000 / 10000 [ 30%] (Warmup)
Chain 4: Iteration: 4000 / 10000 [ 40%] (Warmup)
Chain 4: Iteration: 5000 / 10000 [ 50%] (Warmup)
Chain 4: Iteration: 5001 / 10000 [ 50%] (Sampling)
Chain 4: Iteration: 6000 / 10000 [ 60%] (Sampling)
Chain 4: Iteration: 7000 / 10000 [ 70%] (Sampling)
Chain 4: Iteration: 8000 / 10000 [ 80%] (Sampling)
Chain 4: Iteration: 9000 / 10000 [ 90%] (Sampling)
Chain 4: Iteration: 10000 / 10000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 11.9763 seconds (Warm-up)
Chain 4: 10.1281 seconds (Sampling)
Chain 4: 22.1044 seconds (Total)
Chain 4: mlm_path_plot(dACC.mediation.inverse.results, text = T)
mlm.sum.inverse<-mlm_summary(dACC.mediation.inverse.results, level = .9)
mlm_pars_plot(dACC.mediation.inverse.results, type = "coef", pars = c( "cp", "c", "me","pme"))No id variables; using all as measure variables
# check for convergence of mcmc traces, requires bayesplot package
mcmc_trace(as.data.frame(dACC.mediation.inverse.results), pars = c("a","b","cp","c"))
mcmc_rank_hist(as.data.frame(dACC.mediation.inverse.results), pars = c("a","b","cp","c"))
# Create pretty table of the parameter estimates
formattable(mlm.sum.inverse)| Parameter | Mean | SE | Median | 5% | 95% | n_eff | Rhat |
|---|---|---|---|---|---|---|---|
| a | 0.12 | 0.05 | 0.12 | 0.04 | 0.21 | 15979 | 1 |
| b | 3.45 | 0.77 | 3.45 | 2.20 | 4.71 | 16692 | 1 |
| cp | 2.82 | 0.92 | 2.82 | 1.31 | 4.34 | 8112 | 1 |
| me | 0.45 | 0.26 | 0.44 | 0.06 | 0.89 | 11173 | 1 |
| c | 3.27 | 0.95 | 3.27 | 1.71 | 4.85 | 8351 | 1 |
| pme | 0.14 | 0.24 | 0.13 | 0.02 | 0.30 | 17530 | 1 |
6 Between-Subjects correlations between Reward Rate, Motivation Ratings, and dACC
fit1<-lm(formula = meanRR.normed ~ value.bundle.mean.normed, data = data.dACC.sum.sub)
fit2<-lm(formula = motive_rating.normed ~ value.bundle.mean.normed, data = data.dACC.sum.sub)
fit3<-lm(formula = meanRR.normed ~ value.bundle.mean.normed + motive_rating.normed, data = data.dACC.sum.sub)
tab_model(fit1,fit2,fit3,
title = "Between-Subjects Effects",
dv.labels = c("Reward Rate","Motivation Ratings","Reward Rate"),
pred.labels = c("Intercept","DACC Betas","Motivation Ratings"),
show.re.var = FALSE, show.stat = TRUE, show.icc = FALSE)| Â | Reward Rate | Motivation Ratings | Reward Rate | |||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Predictors | Estimates | CI | Statistic | p | Estimates | CI | Statistic | p | Estimates | CI | Statistic | p |
| Intercept | 0.00 | -0.29 – 0.29 | 0.00 | 1.000 | 0.00 | -0.28 – 0.28 | 0.00 | 1.000 | 0.00 | -0.29 – 0.29 | 0.00 | 1.000 |
| DACC Betas | 0.28 | -0.02 – 0.57 | 1.90 | 0.064 | 0.33 | 0.04 – 0.62 | 2.32 | 0.025 | 0.23 | -0.08 – 0.54 | 1.47 | 0.149 |
| Motivation Ratings | 0.15 | -0.16 – 0.46 | 0.99 | 0.328 | ||||||||
| Observations | 46 | 46 | 46 | |||||||||
| R2 / R2 adjusted | 0.076 / 0.055 | 0.109 / 0.089 | 0.097 / 0.055 | |||||||||
6.1 Manuscript Figure 3b: dACC Activation by Subject (Scatterplot)
6.2 Manuscript Save Figure 3b
p.dACC.bysub<-cowplot::plot_grid(p.dACC.RR.bysub, p.dACC.motive.bysub,
labels = c(" "," "),
rel_widths = c(1, 1),
rel_heights = c(1.15, 1),
align="v", nrow = 1)`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'p.dACC.bysub
# Save the plots
ggsave(filename = "Figure3_dACC_bysubject.png", plot = p.dACC.bysub,
path = figure.path, device = "png", width = 4, height = 2, scale = 2.5, dpi = 300)7 Dorsal ACC Selectively Encodes Subjective Motivational Value and Modulates Motivated Task Performance
7.1 Create table of mixed models of reward rate for forest plot
m.dACC.RR<-lme4::lmer(formula = meanRR ~ mot9Code + value.bundle.mean + (1+mot9Code|subID),
data = data.dACC.sum)
m.vmpfc.RR<-lme4::lmer(formula = meanRR ~ mot9Code + value.bundle.mean + (1+mot9Code|subID),
data = data.vmpfc.sum)
m.accumbens.RR<-lme4::lmer(formula = meanRR ~ mot9Code + value.bundle.mean + (1+mot9Code|subID),
data = data.accumbens.sum)
m.caudate.RR<-lme4::lmer(formula = meanRR ~ mot9Code + value.bundle.mean + (1+mot9Code|subID),
data = data.caudate.sum)
m.putamen.RR<-lme4::lmer(formula = meanRR ~ mot9Code + value.bundle.mean + (1+mot9Code|subID),
data = data.putamen.sum)
m.insula.RR<-lme4::lmer(formula = meanRR ~ mot9Code + value.bundle.mean + (1+mot9Code|subID),
data = data.insula.sum)
# Putting Models in a table
tab_model(m.dACC.RR,m.vmpfc.RR,m.accumbens.RR,m.caudate.RR,m.putamen.RR,m.insula.RR,
title = "Reward Rate Predicted by Bundled Betas",
dv.labels = c("dACC","vmPFC","accumbens","caudate","putamen","anterior insula"),
pred.labels = c("Intercept","Incentive Conditions","Bundled Betas"),
show.re.var = FALSE, show.stat = TRUE, show.icc = FALSE)| Â | dACC | vmPFC | accumbens | caudate | putamen | anterior insula | ||||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Predictors | Estimates | CI | Statistic | p | Estimates | CI | Statistic | p | Estimates | CI | Statistic | p | Estimates | CI | Statistic | p | Estimates | CI | Statistic | p | Estimates | CI | Statistic | p |
| Intercept | 0.68 | 0.65 – 0.71 | 40.93 | <0.001 | 0.70 | 0.66 – 0.73 | 42.70 | <0.001 | 0.69 | 0.66 – 0.73 | 42.99 | <0.001 | 0.69 | 0.66 – 0.72 | 42.93 | <0.001 | 0.69 | 0.66 – 0.72 | 43.19 | <0.001 | 0.69 | 0.66 – 0.73 | 42.87 | <0.001 |
| Incentive Conditions | 0.02 | 0.01 – 0.03 | 2.91 | 0.004 | 0.02 | 0.01 – 0.03 | 3.34 | 0.001 | 0.02 | 0.01 – 0.03 | 3.37 | 0.001 | 0.02 | 0.01 – 0.03 | 3.26 | 0.001 | 0.02 | 0.01 – 0.03 | 3.10 | 0.002 | 0.02 | 0.01 – 0.03 | 3.21 | 0.001 |
| Bundled Betas | 0.09 | 0.02 – 0.17 | 2.41 | 0.016 | 0.01 | -0.03 – 0.06 | 0.55 | 0.585 | -0.00 | -0.02 – 0.01 | -0.58 | 0.562 | 0.03 | -0.02 – 0.08 | 1.06 | 0.288 | 0.06 | -0.01 – 0.13 | 1.68 | 0.093 | 0.02 | -0.06 – 0.10 | 0.51 | 0.607 |
| N | 46 subID | 46 subID | 46 subID | 46 subID | 46 subID | 46 subID | ||||||||||||||||||
| Observations | 414 | 414 | 414 | 414 | 414 | 414 | ||||||||||||||||||
| Marginal R2 / Conditional R2 | 0.040 / 0.591 | 0.026 / 0.598 | 0.026 / 0.596 | 0.028 / 0.593 | 0.031 / 0.594 | 0.026 / 0.595 | ||||||||||||||||||
7.2 Create table of mixed models of motivation ratings for forest plot
m.dACC.motive<-lme4::lmer(formula = motive_rating ~ mot9Code + value.bundle.mean + (1+mot9Code|subID),
data = data.dACC.sum)
m.vmpfc.motive<-lme4::lmer(formula = motive_rating ~ mot9Code + value.bundle.mean + (1+mot9Code|subID),
data = data.vmpfc.sum)
m.accumbens.motive<-lme4::lmer(formula = motive_rating ~ mot9Code + value.bundle.mean + (1+mot9Code|subID),
data = data.accumbens.sum)
m.caudate.motive<-lme4::lmer(formula = motive_rating ~ mot9Code + value.bundle.mean + (1+mot9Code|subID),
data = data.caudate.sum)
m.putamen.motive<-lme4::lmer(formula = motive_rating ~ mot9Code + value.bundle.mean + (1+mot9Code|subID),
data = data.putamen.sum)
m.insula.motive<-lme4::lmer(formula = motive_rating ~ mot9Code + value.bundle.mean + (1+mot9Code|subID),
data = data.insula.sum)
# Putting models into a table
tab_model(m.dACC.motive,m.vmpfc.motive,m.accumbens.motive,m.caudate.motive,m.putamen.motive,m.insula.motive,
title = "Motivation Ratings Predicted by Bundled Betas",
dv.labels = c("dACC","vmPFC","Accumbens","Caudate","Putamen","Anterior Insula"),
pred.labels = c("Intercept","Incentive Conditions","Bundled Betas"),
show.re.var = FALSE, show.stat = TRUE, show.icc = FALSE)| Â | dACC | vmPFC | Accumbens | Caudate | Putamen | Anterior Insula | ||||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Predictors | Estimates | CI | Statistic | p | Estimates | CI | Statistic | p | Estimates | CI | Statistic | p | Estimates | CI | Statistic | p | Estimates | CI | Statistic | p | Estimates | CI | Statistic | p |
| Intercept | 4.70 | 4.45 – 4.94 | 37.43 | <0.001 | 4.90 | 4.67 – 5.13 | 41.04 | <0.001 | 4.93 | 4.70 – 5.16 | 42.06 | <0.001 | 4.87 | 4.64 – 5.11 | 40.65 | <0.001 | 4.87 | 4.64 – 5.09 | 42.32 | <0.001 | 4.87 | 4.64 – 5.10 | 41.89 | <0.001 |
| Incentive Conditions | 0.67 | 0.54 – 0.81 | 9.78 | <0.001 | 0.72 | 0.58 – 0.86 | 10.38 | <0.001 | 0.72 | 0.58 – 0.86 | 10.34 | <0.001 | 0.71 | 0.57 – 0.84 | 10.27 | <0.001 | 0.69 | 0.55 – 0.83 | 9.80 | <0.001 | 0.68 | 0.55 – 0.82 | 9.85 | <0.001 |
| Bundled Betas | 1.64 | 0.83 – 2.45 | 3.96 | <0.001 | -0.32 | -0.84 – 0.21 | -1.19 | 0.233 | -0.02 | -0.19 – 0.15 | -0.22 | 0.825 | 0.74 | 0.11 – 1.37 | 2.29 | 0.022 | 1.29 | 0.51 – 2.07 | 3.22 | 0.001 | 1.36 | 0.43 – 2.29 | 2.87 | 0.004 |
| N | 46 subID | 46 subID | 46 subID | 46 subID | 46 subID | 46 subID | ||||||||||||||||||
| Observations | 414 | 414 | 414 | 414 | 414 | 414 | ||||||||||||||||||
| Marginal R2 / Conditional R2 | 0.294 / 0.501 | 0.261 / 0.494 | 0.258 / 0.491 | 0.267 / 0.499 | 0.279 / 0.502 | 0.275 / 0.497 | ||||||||||||||||||
7.3 Combine Coefficients into Table for Forest Plot
# Reward Rate Coefs
dACC.coef.RR<-coef(summary(m.dACC.RR))
vmpfc.coef.RR<-coef(summary(m.vmpfc.RR))
accumbens.coef.RR<-coef(summary(m.accumbens.RR))
caudate.coef.RR<-coef(summary(m.caudate.RR))
putamen.coef.RR<-coef(summary(m.putamen.RR))
insula.coef.RR<-coef(summary(m.insula.RR))
# Motivation Rating Coefs
dACC.coef.motive<-coef(summary(m.dACC.motive))
vmpfc.coef.motive<-coef(summary(m.vmpfc.motive))
accumbens.coef.motive<-coef(summary(m.accumbens.motive))
caudate.coef.motive<-coef(summary(m.caudate.motive))
putamen.coef.motive<-coef(summary(m.putamen.motive))
insula.coef.motive<-coef(summary(m.insula.motive))
data.SV.names<-rep(c("dACC","vmPFC","NAc","Caudate","Putamen","Insula"),2)
data.SV.betas<-c(dACC.coef.RR["value.bundle.mean","Estimate"],
vmpfc.coef.RR["value.bundle.mean","Estimate"],
accumbens.coef.RR["value.bundle.mean","Estimate"],
caudate.coef.RR["value.bundle.mean","Estimate"],
putamen.coef.RR["value.bundle.mean","Estimate"],
insula.coef.RR["value.bundle.mean","Estimate"],
dACC.coef.motive["value.bundle.mean","Estimate"],
vmpfc.coef.motive["value.bundle.mean","Estimate"],
accumbens.coef.motive["value.bundle.mean","Estimate"],
caudate.coef.motive["value.bundle.mean","Estimate"],
putamen.coef.motive["value.bundle.mean","Estimate"],
insula.coef.motive["value.bundle.mean","Estimate"])
data.SV.std<-c(dACC.coef.RR["value.bundle.mean","Std. Error"],
vmpfc.coef.RR["value.bundle.mean","Std. Error"],
accumbens.coef.RR["value.bundle.mean","Std. Error"],
caudate.coef.RR["value.bundle.mean","Std. Error"],
putamen.coef.RR["value.bundle.mean","Std. Error"],
insula.coef.RR["value.bundle.mean","Std. Error"],
dACC.coef.motive["value.bundle.mean","Std. Error"],
vmpfc.coef.motive["value.bundle.mean","Std. Error"],
accumbens.coef.motive["value.bundle.mean","Std. Error"],
caudate.coef.motive["value.bundle.mean","Std. Error"],
putamen.coef.motive["value.bundle.mean","Std. Error"],
insula.coef.motive["value.bundle.mean","Std. Error"])
data.SV.ci<-data.SV.std*1.96
data.SV.model<-c(rep("Reward Rate",6),rep("Motivation Ratings",6))
data.SV.number<-c(rep(1:6,2))
# Combine into a data frame
data.SV<-data.frame(names=data.SV.names,betas=data.SV.betas,se=data.SV.std,ci=data.SV.ci,Outcome=data.SV.model,number=data.SV.number)
data.SV$names<-factor(data.SV$names, levels = c("vmPFC","NAc","Caudate","Putamen","Insula","dACC"))
# Plot the Bundled Beta Estimates for Each of the ROIs
#p.linecolor<-brewer.pal(n=8, name = 'Set1')
p.bundled<-ggplot(data = data.SV, aes(x = betas, y = names, shape=Outcome, color=Outcome)) +
geom_point(stat="identity", size = 3) +
geom_errorbarh(aes(xmin=betas-ci, xmax=betas+ci, height = .1)) +
geom_vline(xintercept=0, linetype="dashed", color="grey0", size = .7) +
xlab("Bundled Beta Effect") + ylab("Regions of Interest") +
#xlim(-0.05, .15) +
facet_grid(.~Outcome, scales = "free_x")+
scale_color_brewer(palette="Set1") +
theme_classic() +
theme(#plot.title=element_text(size=22,face="bold", vjust=2),
axis.title=element_text(size=14, face = "bold"),
axis.text=element_text(size=14),
legend.position="none",
strip.text.x = element_text(size = 14, face = "bold"))
p.bundled
7.4 Manuscript Save Figure 4: Forest Plot
p.figure4<-cowplot::plot_grid(p.bundled,
labels = c("4"), label_size = 16)
ggsave(filename = "Figure4.eps", plot = p.figure4,
device="eps", dpi = 300,
path = figure.path, width = 7, height = 5, scale = 1.1)8 Exploratory Visulization and Plots (Not included in publication)
8.1 ROI Analyses: vmPFC
p.vmpfc<-cowplot::plot_grid(p.vmpfc.RR.bysub, p.vmpfc.motive.bysub,
labels = c(" "," "),
rel_widths = c(1, 1),
rel_heights = c(1.15, 1),
align="v", nrow = 1)`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'p.vmpfc
8.1.1 Within-Subjects Scatterplots vmPFC
p.vmpfc.diff<-cowplot::plot_grid(p.vmpfc.MotDiff, p.vmpfc.MotDiff2,
labels = c(" "," "),
rel_widths = c(1, 1),
rel_heights = c(1, 1),
align="v", nrow = 1)`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'p.vmpfc.diff
8.2 ROI Analyses: Accumbens
p.accumbens<-cowplot::plot_grid(p.accumbens.RR.bysub, p.accumbens.motive.bysub,
labels = c(" "," "),
rel_widths = c(1, 1),
rel_heights = c(1, 1),
align="v", nrow = 1)`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'p.accumbens
8.3 ROI Analyses: Caudate
p.caudate<-cowplot::plot_grid(p.caudate.RR.bysub, p.caudate.motive.bysub,
labels = c(" "," "),
rel_widths = c(1, 1),
rel_heights = c(1.15, 1),
align="v", nrow = 1)`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'p.caudate
8.4 Within-Subjects Scatterplots Caudate
p.caudate.diff<-cowplot::plot_grid(p.caudate.MotDiff, p.caudate.MotDiff2,
labels = c(" "," "),
rel_widths = c(1, 1),
rel_heights = c(1, 1),
align="v", nrow = 1)`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'p.caudate.diff
9 ROI Analyses: Putamen
p.putamen<-cowplot::plot_grid(p.putamen.RR.bysub, p.putamen.motive.bysub,
labels = c(" "," "),
rel_widths = c(1, 1),
rel_heights = c(1, 1),
align="v", nrow = 1)`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'p.putamen
9.1 Within-Subjects Scatterplots Putamen
p.putamen.diff<-cowplot::plot_grid(p.putamen.MotDiff, p.putamen.MotDiff2,
labels = c(" "," "),
rel_widths = c(1, 1),
rel_heights = c(1, 1),
align="v", nrow = 1)`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'p.putamen.diff
9.2 Exploratory Visualization: Correlation Plot of bundled betas (by averaged across all conditions within a subject)
data.corr <- data.dACC.sum.sub %>% dplyr::select(subID, dACC = value.bundle.mean) %>%
left_join(data.vmpfc.sum.sub, by = "subID") %>% dplyr::select(subID, dACC, vmpfc = value.bundle.mean) %>%
left_join(data.accumbens.sum.sub, by = "subID") %>% dplyr::select(subID, dACC, vmpfc, accumbens = value.bundle.mean) %>%
left_join(data.caudate.sum.sub, by = "subID") %>% dplyr::select(subID, dACC, vmpfc, accumbens, caudate = value.bundle.mean) %>%
left_join(data.putamen.sum.sub, by = "subID") %>% dplyr::select(subID, dACC, vmpfc, accumbens, caudate, putamen = value.bundle.mean) %>%
left_join(data.insula.sum.sub, by = "subID") %>% dplyr::select(subID, dACC, vmpfc, accumbens, caudate, putamen,insula = value.bundle.mean) %>%
dplyr::select(-subID)
ggcorrplot(corr = cor(data.corr), hc.order=TRUE, type = "lower",
colors = c("#6D9EC1", "white", "#E46726"),
lab = TRUE,
title = "Correlation Matrix of Bundled Beta Estimates Across ROIs",
outline.color = "black", show.diag = TRUE,
p.mat = cor_pmat(data.corr))
10 Exploratory Visualization: Congruency
RR.sum=summarySEwithin2(data=incentive.means.cong, measurevar = "meanRR",
withinvars = c("congruency"), idvar = "subID")Automatically converting the following non-factors to factors: congruencyp.RR.cong<-ggplot(RR.sum, aes(x=congruency, y=meanRR)) +
geom_bar(position=position_dodge(width=0.8), color="black",
stat="identity", width=0.8) +
geom_errorbar(position=position_dodge(width=0.8),
aes(ymin=meanRR-se, ymax=meanRR+se), width=.2) +
xlab("Congruency") + ylab("Reward Rate") +
labs(fill = "Liquid") +
geom_hline(yintercept = .4, linetype="dashed") +
scale_fill_brewer(palette="Set2") +
coord_cartesian(ylim=c(.5,.75)) +
#coord_cartesian(ylim=c(600,650)) +
#scale_fill_discrete(name="Monetary Reward") +
theme(#plot.title=element_text(size=22,face="bold", vjust=2),
axis.title=element_text(size=14,face = "bold"),
axis.text=element_text(size=14),
legend.position="right",
legend.title=element_text(size=14, face = "bold"),
legend.text = element_text(size=14),
strip.text.x = element_text(size = 12))
p.RR.cong
RR.sum=summarySEwithin2(data=incentive9.means.cong, measurevar = "meanRR",
withinvars = c("congruency","Liquid","MoneyLabel"), idvar = "subID")Automatically converting the following non-factors to factors: congruencyp.RR.mot9.cong<-ggplot(RR.sum, aes(x=MoneyLabel, y=meanRR, fill=Liquid)) +
geom_bar(position=position_dodge(width=0.8), color="black",
stat="identity", width=0.8) +
geom_errorbar(position=position_dodge(width=0.8),
aes(ymin=meanRR-se, ymax=meanRR+se), width=.2) +
xlab("Money") + ylab("Reward Rate") +
labs(fill = "Liquid") +
geom_hline(yintercept = .4, linetype="dashed") +
scale_fill_brewer(palette="Set2") +
coord_cartesian(ylim=c(.5,.8)) +
facet_grid(~congruency)+
#scale_fill_discrete(name="Monetary Reward") +
theme(#plot.title=element_text(size=22,face="bold", vjust=2),
axis.title=element_text(size=14,face = "bold"),
axis.text=element_text(size=14),
legend.position="right",
legend.title=element_text(size=14, face = "bold"),
legend.text = element_text(size=14),
strip.text.x = element_text(size = 12))
p.RR.mot9.cong
ggsave(filename = "congruency.eps", plot = p.RR.cong,
device="eps", dpi = 300,
path = figure.path, width = 3, height = 3, scale = 1.1)
ggsave(filename = "congruency_mot9.eps", plot = p.RR.mot9.cong,
device="eps", dpi = 300,
path = figure.path, width = 6, height = 3, scale = 1.1)11 Exploratory Visualization: Reward Rate switch costs
RR.sum=summarySEwithin2(data=incentive.means.switch, measurevar = "meanRR",
withinvars = c("switch_trial"), idvar = "subID")Automatically converting the following non-factors to factors: switch_trialp.RR.switch<-ggplot(RR.sum, aes(x=switch_trial, y=meanRR)) +
geom_bar(position=position_dodge(width=0.8), color="black",
stat="identity", width=0.8) +
geom_errorbar(position=position_dodge(width=0.8),
aes(ymin=meanRR-se, ymax=meanRR+se), width=.2) +
#xlab("Congruency") +
ylab("Reward Rate") +
labs(fill = "Liquid") +
geom_hline(yintercept = .4, linetype="dashed") +
scale_fill_brewer(palette="Set2") +
coord_cartesian(ylim=c(.5,.75)) +
#coord_cartesian(ylim=c(600,650)) +
#scale_fill_discrete(name="Monetary Reward") +
theme(#plot.title=element_text(size=22,face="bold", vjust=2),
axis.title=element_text(size=14,face = "bold"),
axis.text=element_text(size=14),
legend.position="right",
legend.title=element_text(size=14, face = "bold"),
legend.text = element_text(size=14),
strip.text.x = element_text(size = 12))
p.RR.switch
RR.sum=summarySEwithin2(data=incentive9.means.switch, measurevar = "meanRR",
withinvars = c("switch_trial","Liquid","MoneyLabel"), idvar = "subID")Automatically converting the following non-factors to factors: switch_trialp.RR.mot9.switch<-ggplot(RR.sum, aes(x=MoneyLabel, y=meanRR, fill=Liquid)) +
geom_bar(position=position_dodge(width=0.8), color="black",
stat="identity", width=0.8) +
geom_errorbar(position=position_dodge(width=0.8),
aes(ymin=meanRR-se, ymax=meanRR+se), width=.2) +
xlab("Money") + ylab("Reward Rate") +
labs(fill = "Liquid") +
geom_hline(yintercept = .4, linetype="dashed") +
scale_fill_brewer(palette="Set2") +
coord_cartesian(ylim=c(.5,.8)) +
facet_grid(~switch_trial)+
#scale_fill_discrete(name="Monetary Reward") +
theme(#plot.title=element_text(size=22,face="bold", vjust=2),
axis.title=element_text(size=14,face = "bold"),
axis.text=element_text(size=14),
legend.position="right",
legend.title=element_text(size=14, face = "bold"),
legend.text = element_text(size=14),
strip.text.x = element_text(size = 12))
p.RR.mot9.switch
ggsave(filename = "switchcostRR.eps", plot = p.RR.switch,
device="eps", dpi = 300,
path = figure.path, width = 3, height = 3, scale = 1.1)
ggsave(filename = "switchcostRR_mot9.eps", plot = p.RR.mot9.switch,
device="eps", dpi = 300,
path = figure.path, width = 6, height = 3, scale = 1.1)12 Exploratory Visualization: RT Switch Costs
RT.sum=summarySEwithin2(data=incentiveRT.means.switch, measurevar = "meanRT",
withinvars = c("switch_trial"), idvar = "subID")Automatically converting the following non-factors to factors: switch_trialp.RT.switch<-ggplot(RT.sum, aes(x=switch_trial, y=meanRT)) +
geom_bar(position=position_dodge(width=0.8), color="black",
stat="identity", width=0.8) +
geom_errorbar(position=position_dodge(width=0.8),
aes(ymin=meanRT-se, ymax=meanRT+se), width=.2) +
#xlab("Congruency") +
ylab("Response Time") +
labs(fill = "Liquid") +
geom_hline(yintercept = .4, linetype="dashed") +
scale_fill_brewer(palette="Set2") +
#coord_cartesian(ylim=c(.5,.75)) +
coord_cartesian(ylim=c(600,700)) +
#scale_fill_discrete(name="Monetary Reward") +
theme(#plot.title=element_text(size=22,face="bold", vjust=2),
axis.title=element_text(size=14,face = "bold"),
axis.text=element_text(size=14),
legend.position="right",
legend.title=element_text(size=14, face = "bold"),
legend.text = element_text(size=14),
strip.text.x = element_text(size = 12))
p.RT.switch
RT.sum=summarySEwithin2(data=incentive9RT.means.switch, measurevar = "meanRT",
withinvars = c("switch_trial","Liquid","MoneyLabel"), idvar = "subID")Automatically converting the following non-factors to factors: switch_trialp.RT.mot9.switch<-ggplot(RT.sum, aes(x=MoneyLabel, y=meanRT, fill=Liquid)) +
geom_bar(position=position_dodge(width=0.8), color="black",
stat="identity", width=0.8) +
geom_errorbar(position=position_dodge(width=0.8),
aes(ymin=meanRT-se, ymax=meanRT+se), width=.2) +
xlab("Money") + ylab("Response Time") +
labs(fill = "Liquid") +
geom_hline(yintercept = .4, linetype="dashed") +
scale_fill_brewer(palette="Set2") +
#coord_cartesian(ylim=c(.5,.8)) +
coord_cartesian(ylim=c(550,700)) +
facet_grid(~switch_trial)+
#scale_fill_discrete(name="Monetary Reward") +
theme(#plot.title=element_text(size=22,face="bold", vjust=2),
axis.title=element_text(size=14,face = "bold"),
axis.text=element_text(size=14),
legend.position="right",
legend.title=element_text(size=14, face = "bold"),
legend.text = element_text(size=14),
strip.text.x = element_text(size = 12))
p.RT.mot9.switch
#m.switch<-lmer(meanRT ~ switch_trial * liqCode * moneyCode + (1|subID), data = incentive9RT.means.switch)
#summary(m.switch)
ggsave(filename = "switchcostRT.eps", plot = p.RT.switch,
device="eps", dpi = 300,
path = figure.path, width = 3, height = 3, scale = 1.1)
ggsave(filename = "switchcostRT_mot9.eps", plot = p.RT.mot9.switch,
device="eps", dpi = 300,
path = figure.path, width = 6, height = 3, scale = 1.1)